|
0
|
1 /*
|
|
|
2 * Copyright 1996-2006 Catherine Loader.
|
|
|
3 */
|
|
|
4
|
|
|
5 #include "mex.h"
|
|
|
6 /*
|
|
|
7 * Copyright 1996-2006 Catherine Loader.
|
|
|
8 */
|
|
|
9 /*
|
|
|
10 * Copyright (c) 1998-2006 Catherine Loader.
|
|
|
11 * This file contains functions to compute the constants
|
|
|
12 * appearing in the tube formula.
|
|
|
13 */
|
|
|
14
|
|
|
15 #include <stdio.h>
|
|
|
16 #include <math.h>
|
|
|
17 #include "tube.h"
|
|
|
18
|
|
|
19 static double *fd, *ft;
|
|
|
20 static int globm, (*wdf)(), use_covar, kap_terms;
|
|
|
21
|
|
|
22 int k0_reqd(d,n,uc)
|
|
|
23 int d, n, uc;
|
|
|
24 { int m;
|
|
|
25 m = d*(d+1)+1;
|
|
|
26 if (uc) return(2*m*m);
|
|
|
27 else return(2*n*m);
|
|
|
28 }
|
|
|
29
|
|
|
30 void assignk0(z,d,n) /* z should be n*(2*d*d+2*d+2); */
|
|
|
31 double *z;
|
|
|
32 int d, n;
|
|
|
33 { ft = z; z += n*(d*(d+1)+1);
|
|
|
34 fd = z; z += n*(d*(d+1)+1);
|
|
|
35 }
|
|
|
36
|
|
|
37 /* Residual projection of y to the columns of A,
|
|
|
38 * (I - A(R^TR)^{-1}A^T)y
|
|
|
39 * R should be from the QR-decomp. of A.
|
|
|
40 */
|
|
|
41 void rproject(y,A,R,n,p)
|
|
|
42 double *y, *A, *R;
|
|
|
43 int n, p;
|
|
|
44 { double v[1+TUBE_MXDIM];
|
|
|
45 int i, j;
|
|
|
46
|
|
|
47 for (i=0; i<p; i++) v[i] = innerprod(&A[i*n],y,n);
|
|
|
48 qrsolv(R,v,n,p);
|
|
|
49 for (i=0; i<n; i++)
|
|
|
50 for (j=0; j<p; j++)
|
|
|
51 y[i] -= A[j*n+i]*v[j];
|
|
|
52 }
|
|
|
53
|
|
|
54 double k2c(lij,A,m,dd,d)
|
|
|
55 double *lij, *A;
|
|
|
56 int m, d, dd;
|
|
|
57 { int i, j, k, l;
|
|
|
58 double sum, *bk, v[TUBE_MXDIM];
|
|
|
59
|
|
|
60 for (i=0; i<dd*d; i++)
|
|
|
61 chol_hsolve(fd,&lij[i*m],m,dd+1);
|
|
|
62 for (i=0; i<dd*d; i++)
|
|
|
63 for (j=0; j<dd*d; j++)
|
|
|
64 lij[i*m+j+d+1] -= innerprod(&lij[i*m],&lij[j*m],dd+1);
|
|
|
65
|
|
|
66 sum = 0;
|
|
|
67 for (i=0; i<dd; i++)
|
|
|
68 for (j=0; j<i; j++)
|
|
|
69 { bk = &lij[i*d*m + j*d + d+1];
|
|
|
70 for (k=0; k<dd; k++)
|
|
|
71 { v[0] = 0;
|
|
|
72 for (l=0; l<dd; l++) v[l+1] = bk[k*m+l];
|
|
|
73 chol_solve(fd,v,m,dd+1);
|
|
|
74 for (l=0; l<dd; l++) bk[k*m+l] = v[l+1];
|
|
|
75 }
|
|
|
76 for (k=0; k<dd; k++)
|
|
|
77 { v[0] = 0;
|
|
|
78 for (l=0; l<dd; l++) v[l+1] = bk[l*m+k];
|
|
|
79 chol_solve(fd,v,m,dd+1);
|
|
|
80 for (l=0; l<dd; l++) bk[l*m+k] = v[l+1];
|
|
|
81 }
|
|
|
82 sum += bk[i*m+j] - bk[j*m+i];
|
|
|
83 }
|
|
|
84 return(sum*fd[0]*fd[0]);
|
|
|
85 }
|
|
|
86
|
|
|
87 double k2x(lij,A,m,d,dd)
|
|
|
88 double *lij, *A;
|
|
|
89 int m, d, dd;
|
|
|
90 { int i, j, k;
|
|
|
91 double s, v[1+TUBE_MXDIM], *ll;
|
|
|
92
|
|
|
93 /* residual projections of lij onto A = [l,l1,...,ld] */
|
|
|
94 for (i=0; i<d; i++)
|
|
|
95 for (j=i; j<d; j++)
|
|
|
96 { ll = &lij[(i*dd+j)*m];
|
|
|
97 rproject(ll,A,fd,m,d+1);
|
|
|
98 if (i!=j) memcpy(&lij[(j*dd+i)*m],ll,m*sizeof(double));
|
|
|
99 }
|
|
|
100
|
|
|
101 /* compute lij[j][i] = e_i^T (A^T A)^{-1} B_j^T */
|
|
|
102 for (k=0; k<m; k++)
|
|
|
103 for (j=0; j<d; j++)
|
|
|
104 { v[0] = 0;
|
|
|
105 for (i=0; i<d; i++) v[i+1] = lij[(j*dd+i)*m+k];
|
|
|
106 qrsolv(fd,v,m,d+1);
|
|
|
107 for (i=0; i<d; i++) lij[(j*dd+i)*m+k] = v[i+1];
|
|
|
108 }
|
|
|
109
|
|
|
110 /* finally, add up to get the kappa2 term */
|
|
|
111 s = 0;
|
|
|
112 for (j=0; j<d; j++)
|
|
|
113 for (k=0; k<j; k++)
|
|
|
114 s += innerprod(&lij[(j*dd+j)*m],&lij[(k*dd+k)*m],m)
|
|
|
115 - innerprod(&lij[(j*dd+k)*m],&lij[(k*dd+j)*m],m);
|
|
|
116
|
|
|
117 return(s*fd[0]*fd[0]);
|
|
|
118 }
|
|
|
119
|
|
|
120 void d2c(ll,nn,li,ni,lij,nij,M,m,dd,d)
|
|
|
121 double *ll, *nn, *li, *ni, *lij, *nij, *M;
|
|
|
122 int m, dd, d;
|
|
|
123 { int i, j, k, l, t, u, v, w;
|
|
|
124 double z;
|
|
|
125
|
|
|
126 for (i=0; i<dd; i++)
|
|
|
127 for (j=0; j<dd; j++)
|
|
|
128 { for (k=0; k<d; k++)
|
|
|
129 { for (l=0; l<d; l++)
|
|
|
130 { z = M[i*d+k]*M[j*d+l];
|
|
|
131 if (z != 0.0)
|
|
|
132 { nij[(i*d+j)*m] += z*lij[(k*d+l)*m];
|
|
|
133 for (t=0; t<d; t++) /* need d, not dd here */
|
|
|
134 for (u=0; u<d; u++)
|
|
|
135 nij[(i*d+j)*m+t+1] += z*M[t*d+u]*lij[(k*d+l)*m+u+1];
|
|
|
136 for (t=0; t<dd; t++)
|
|
|
137 for (u=0; u<dd; u++)
|
|
|
138 { for (v=0; v<d; v++)
|
|
|
139 for (w=0; w<d; w++)
|
|
|
140 nij[(i*d+j)*m+(t*d+u)+d+1] +=
|
|
|
141 z*M[t*d+v]*M[u*d+w]*lij[(k*d+l)*m+(v*d+w)+d+1];
|
|
|
142 for (v=0; v<d; v++)
|
|
|
143 nij[(i*d+j)*m+(t*d+u)+d+1] += z*M[(v+1)*d*d+t*d+u]*lij[(k*d+l)*m+v+1];
|
|
|
144 }
|
|
|
145 }
|
|
|
146 }
|
|
|
147
|
|
|
148 z = M[(k+1)*d*d+i*d+j];
|
|
|
149 if (z!=0.0)
|
|
|
150 { nij[(i*d+j)*m] += z*li[k*m];
|
|
|
151 for (t=0; t<d; t++)
|
|
|
152 for (u=0; u<d; u++)
|
|
|
153 nij[(i*d+j)*m+t+1] += z*M[t*d+u]*li[k*m+u+1];
|
|
|
154 for (t=0; t<dd; t++)
|
|
|
155 for (u=0; u<dd; u++)
|
|
|
156 { for (v=0; v<d; v++)
|
|
|
157 for (w=0; w<d; w++)
|
|
|
158 nij[(i*d+j)*m+(t*d+u)+d+1] += z*M[t*d+v]*M[u*d+w]*lij[(v*d+w)*m+k+1];
|
|
|
159 for (v=0; v<d; v++)
|
|
|
160 nij[(i*d+j)*m+(t*d+u)+d+1] += z*M[(v+1)*d*d+t*d+u]*li[k*m+v+1];
|
|
|
161 }
|
|
|
162 }
|
|
|
163 }
|
|
|
164 }
|
|
|
165 }
|
|
|
166
|
|
|
167 void d2x(li,lij,nij,M,m,dd,d)
|
|
|
168 double *li, *lij, *nij, *M;
|
|
|
169 int m, dd, d;
|
|
|
170 { int i, j, k, l, z;
|
|
|
171 double t;
|
|
|
172 for (i=0; i<dd; i++)
|
|
|
173 for (j=0; j<dd; j++)
|
|
|
174 { for (k=0; k<d; k++)
|
|
|
175 { for (l=0; l<d; l++)
|
|
|
176 { t = M[i*d+k] * M[j*d+l];
|
|
|
177 if (t != 0.0)
|
|
|
178 { for (z=0; z<m; z++)
|
|
|
179 nij[(i*d+j)*m+z] += t*lij[(k*d+l)*m+z];
|
|
|
180 }
|
|
|
181 }
|
|
|
182 t = M[(k+1)*d*d+i*d+j];
|
|
|
183 if (t!=0.0)
|
|
|
184 for (z=0; z<m; z++)
|
|
|
185 nij[(i*d+j)*m+z] += t*li[k*m+z];
|
|
|
186 }
|
|
|
187 }
|
|
|
188 }
|
|
|
189
|
|
|
190 int k0x(x,d,kap,M)
|
|
|
191 double *x, *kap, *M;
|
|
|
192 int d;
|
|
|
193 { double det, *lij, *nij, z;
|
|
|
194 int j, m, r;
|
|
|
195
|
|
|
196 r = 1 + ((d>=2) & (kap_terms >= 3));
|
|
|
197 m = globm = wdf(x,ft,r);
|
|
|
198
|
|
|
199 memcpy(fd,ft,m*(d+1)*sizeof(double));
|
|
|
200 if (use_covar) chol_dec(fd,m,d+1);
|
|
|
201 else qr(fd,m,d+1,NULL);
|
|
|
202
|
|
|
203 det = 1;
|
|
|
204 for (j=1; j<=d; j++)
|
|
|
205 det *= fd[j*(m+1)]/fd[0];
|
|
|
206 kap[0] = det;
|
|
|
207 if (kap_terms == 1) return(1);
|
|
|
208 kap[1] = 0.0;
|
|
|
209 if ((kap_terms == 2) | (d<=1)) return(2);
|
|
|
210
|
|
|
211 lij = &ft[(d+1)*m];
|
|
|
212 nij = &fd[(d+1)*m];
|
|
|
213 memcpy(nij,lij,m*d*d*sizeof(double));
|
|
|
214 z = (use_covar) ? k2c(nij,ft,m,d,d) : k2x(nij,ft,m,d,d);
|
|
|
215 kap[2] = z*det;
|
|
|
216 if ((kap_terms == 3) | (d==2)) return(3);
|
|
|
217
|
|
|
218 kap[3] = 0;
|
|
|
219 return(4);
|
|
|
220 }
|
|
|
221
|
|
|
222 void d1c(li,ni,m,d,M)
|
|
|
223 double *li, *ni, *M;
|
|
|
224 int m, d;
|
|
|
225 { int i, j, k, l;
|
|
|
226 double t;
|
|
|
227
|
|
|
228 fd[0] = ft[0];
|
|
|
229 for (i=0; i<d; i++)
|
|
|
230 { t = 0;
|
|
|
231 for (j=0; j<d; j++) t += M[i*d+j]*li[j*m];
|
|
|
232 fd[i+1] = ni[i*m] = t;
|
|
|
233
|
|
|
234 for (j=0; j<d; j++)
|
|
|
235 { t = 0;
|
|
|
236 for (k=0; k<d; k++)
|
|
|
237 for (l=0; l<d; l++)
|
|
|
238 t += li[k*m+l+1] * M[i*d+k] * M[j*d+l];
|
|
|
239 ni[i*m+j+1] = t;
|
|
|
240 }
|
|
|
241 }
|
|
|
242 }
|
|
|
243
|
|
|
244 void d1x(li,ni,m,d,M)
|
|
|
245 double *li, *ni, *M;
|
|
|
246 int m, d;
|
|
|
247 { int i, j, k;
|
|
|
248 memcpy(fd,ft,m*sizeof(double));
|
|
|
249 setzero(ni,m*d);
|
|
|
250 for (j=0; j<d; j++)
|
|
|
251 for (k=0; k<d; k++)
|
|
|
252 if (M[j*d+k]!=0)
|
|
|
253 for (i=0; i<m; i++) ni[j*m+i] += M[j*d+k]*li[k*m+i];
|
|
|
254 }
|
|
|
255
|
|
|
256 int l1x(x,d,lap,M)
|
|
|
257 double *x, *lap, *M;
|
|
|
258 int d;
|
|
|
259 { double det, sumcj, *u, v[TUBE_MXDIM];
|
|
|
260 double *ll, *li, *lij, *ni, *nij;
|
|
|
261 int i, j, m;
|
|
|
262 if (kap_terms<=1) return(0);
|
|
|
263 m = globm;
|
|
|
264 li = &ft[m]; lij = &ft[(d+1)*m];
|
|
|
265 ni = &fd[m]; nij = &fd[(d+1)*m];
|
|
|
266 setzero(ni,m*d);
|
|
|
267 setzero(nij,m*d*d);
|
|
|
268
|
|
|
269 if (use_covar) d1c(li,ni,m,d,M);
|
|
|
270 else d1x(li,ni,m,d,M);
|
|
|
271
|
|
|
272 /* the last (d+1) columns of nij are free, use for an extra copy of ni */
|
|
|
273 ll = &fd[d*d*m];
|
|
|
274 u = &ll[d*m];
|
|
|
275 if (use_covar)
|
|
|
276 memcpy(u,&ni[(d-1)*m],d*sizeof(double)); /* cov(ld, (l,l1,...ld-1)) */
|
|
|
277 else
|
|
|
278 memcpy(ll,fd,(d+1)*m*sizeof(double));
|
|
|
279
|
|
|
280 if (use_covar) chol_dec(fd,m,d+1);
|
|
|
281 else qr(fd,m,d+1,NULL);
|
|
|
282 det = 1;
|
|
|
283 for (j=1; j<d; j++)
|
|
|
284 det *= fd[(m+1)*j]/fd[0];
|
|
|
285 lap[0] = det;
|
|
|
286 if ((kap_terms==2) | (d<=1)) return(1);
|
|
|
287
|
|
|
288 sumcj = 0.0;
|
|
|
289 if (use_covar)
|
|
|
290 { d2c(ft,fd,li,ni,lij,nij,M,m,d-1,d);
|
|
|
291 chol_solve(fd,u,m,d);
|
|
|
292 for (i=0; i<d-1; i++)
|
|
|
293 { v[0] = 0;
|
|
|
294 for (j=0; j<d-1; j++)
|
|
|
295 v[j+1] = nij[(i*d+j)*m+d] - innerprod(u,&nij[(i*d+j)*m],d);
|
|
|
296 chol_solve(fd,v,m,d);
|
|
|
297 sumcj -= v[i+1];
|
|
|
298 }
|
|
|
299 }
|
|
|
300 else
|
|
|
301 { d2x(li,lij,nij,M,m,d-1,d);
|
|
|
302 rproject(u,ll,fd,m,d);
|
|
|
303 for (i=0; i<d-1; i++)
|
|
|
304 { v[0] = 0;
|
|
|
305 for (j=0; j<d-1; j++) v[j+1] = innerprod(&nij[(i*d+j)*m],u,m);
|
|
|
306 qrsolv(fd,v,m,d);
|
|
|
307 sumcj -= v[i+1];
|
|
|
308 }
|
|
|
309 }
|
|
|
310
|
|
|
311 lap[1] = sumcj*det*fd[0]/fd[(m+1)*d];
|
|
|
312 if ((kap_terms==3) | (d==2)) return(2);
|
|
|
313
|
|
|
314 if (use_covar) lap[2] = k2c(nij,ll,m,d-1,d)*det;
|
|
|
315 else lap[2] = k2x(nij,ll,m,d-1,d)*det;
|
|
|
316 return(3);
|
|
|
317 }
|
|
|
318
|
|
|
319 int m0x(x,d,m0,M)
|
|
|
320 double *x, *m0, *M;
|
|
|
321 int d;
|
|
|
322 { double det, *li, *ni, *lij, *nij, *ll, *u1, *u2;
|
|
|
323 double om, so, co, sumcj, v[TUBE_MXDIM];
|
|
|
324 int m, i, j;
|
|
|
325
|
|
|
326 if ((kap_terms<=2) | (d<=1)) return(0);
|
|
|
327
|
|
|
328 m = globm;
|
|
|
329 li = &ft[m]; lij = &ft[(d+1)*m];
|
|
|
330 ni = &fd[m]; nij = &fd[(d+1)*m];
|
|
|
331 setzero(ni,m*d);
|
|
|
332 setzero(nij,m*d*d);
|
|
|
333
|
|
|
334 if (use_covar) d1c(li,ni,m,d,M);
|
|
|
335 else d1x(li,ni,m,d,M);
|
|
|
336
|
|
|
337 /* the last (d+1) columns of nij are free, use for an extra copy of ni */
|
|
|
338 ll = &fd[d*d*m];
|
|
|
339 u1 = &ll[d*m];
|
|
|
340 u2 = &ll[(d-1)*m];
|
|
|
341 if (use_covar)
|
|
|
342 { memcpy(u1,&ni[(d-1)*m],d*sizeof(double));
|
|
|
343 memcpy(u2,&ni[(d-2)*m],d*sizeof(double));
|
|
|
344 }
|
|
|
345 else
|
|
|
346 memcpy(ll,fd,(d+1)*m*sizeof(double));
|
|
|
347
|
|
|
348 if (use_covar) chol_dec(fd,m,d+1);
|
|
|
349 else qr(fd,m,d+1,NULL);
|
|
|
350 det = 1;
|
|
|
351 for (j=1; j<d-1; j++)
|
|
|
352 det *= fd[j*(m+1)]/fd[0];
|
|
|
353 om = atan2(fd[d*(m+1)],-fd[d*(m+1)-1]);
|
|
|
354 m0[0] = det*om;
|
|
|
355 if ((kap_terms==3) | (d==2)) return(1);
|
|
|
356
|
|
|
357 so = sin(om)/fd[d*(m+1)];
|
|
|
358 co = (1-cos(om))/fd[(d-1)*(m+1)];
|
|
|
359 sumcj = 0.0;
|
|
|
360 if (use_covar)
|
|
|
361 { d2c(ft,fd,li,ni,lij,nij,M,m,d-2,d);
|
|
|
362 chol_solve(fd,u1,m,d);
|
|
|
363 chol_solve(fd,u2,m,d-1);
|
|
|
364 for (i=0; i<d-2; i++)
|
|
|
365 { v[0] = 0;
|
|
|
366 for (j=0; j<d-2; j++)
|
|
|
367 v[j+1] =
|
|
|
368 so*(nij[(i*d+j)*m+d]-innerprod(u1,&nij[(i*d+j)*m],d))
|
|
|
369 +co*(nij[(i*d+j)*m+d-1]-innerprod(u2,&nij[(i*d+j)*m],d-1));
|
|
|
370 qrsolv(fd,v,m,d-1);
|
|
|
371 sumcj -= v[i+1];
|
|
|
372 }
|
|
|
373 }
|
|
|
374 else
|
|
|
375 { d2x(li,lij,nij,M,m,d-2,d);
|
|
|
376 rproject(u1,ll,fd,m,d);
|
|
|
377 rproject(u2,ll,fd,m,d-1); /* now, u1, u2 are unnormalized n1*, n2* */
|
|
|
378 for (i=0; i<m; i++)
|
|
|
379 u1[i] = so*u1[i] + co*u2[i]; /* for n1*, n2* */
|
|
|
380 for (i=0; i<d-2; i++)
|
|
|
381 { v[0] = 0;
|
|
|
382 for (j=0; j<d-2; j++)
|
|
|
383 v[j+1] = innerprod(&nij[(i*d+j)*m],u1,m);
|
|
|
384 qrsolv(fd,v,m,d-1);
|
|
|
385 sumcj -= v[i+1];
|
|
|
386 }
|
|
|
387 }
|
|
|
388
|
|
|
389 m0[1] = sumcj*det*fd[0];
|
|
|
390 return(2);
|
|
|
391 }
|
|
|
392
|
|
|
393 int n0x(x,d,n0,M)
|
|
|
394 double *x, *n0, *M;
|
|
|
395 int d;
|
|
|
396 { double det, *li, *ni, *a0, *a1, *a2;
|
|
|
397 int j, m;
|
|
|
398
|
|
|
399 if ((kap_terms <= 3) | (d <= 2)) return(0);
|
|
|
400
|
|
|
401 m = globm;
|
|
|
402 li = &ft[m];
|
|
|
403 ni = &fd[m];
|
|
|
404
|
|
|
405 if (use_covar) d1c(li,ni,m,d,M);
|
|
|
406 else d1x(li,ni,m,d,M);
|
|
|
407
|
|
|
408 det = 1;
|
|
|
409 if (use_covar) chol_dec(fd,m,d+1);
|
|
|
410 else qr(fd,m,d+1,NULL);
|
|
|
411 for (j=1; j<d-2; j++)
|
|
|
412 det *= fd[j*(m+1)]/fd[0];
|
|
|
413
|
|
|
414 a0 = &ni[(d-3)*m+d-2];
|
|
|
415 a1 = &ni[(d-2)*m+d-2];
|
|
|
416 a2 = &ni[(d-1)*m+d-2];
|
|
|
417
|
|
|
418 a0[0] = a1[1]*a2[2];
|
|
|
419 a0[1] =-a1[0]*a2[2];
|
|
|
420 a0[2] = a1[0]*a2[1]-a1[1]*a2[0];
|
|
|
421 a1[0] = 0;
|
|
|
422 a1[1] = a2[2];
|
|
|
423 a1[2] =-a2[1];
|
|
|
424 a2[0] = a2[1] = 0.0; a2[2] = 1.0;
|
|
|
425 rn3(a0); rn3(a1);
|
|
|
426 n0[0] = det*sptarea(a0,a1,a2);
|
|
|
427 return(1);
|
|
|
428 }
|
|
|
429
|
|
|
430 int kodf(ll,ur,mg,kap,lap)
|
|
|
431 double *ll, *ur, *kap, *lap;
|
|
|
432 int *mg;
|
|
|
433 { double x[1], *l0, *l1, t, sum;
|
|
|
434 int i, j, n;
|
|
|
435
|
|
|
436 sum = 0.0;
|
|
|
437 for (i=0; i<=mg[0]; i++)
|
|
|
438 { if (i&1) { l1 = fd; l0 = ft; }
|
|
|
439 else { l1 = ft; l0 = fd; }
|
|
|
440 x[0] = ll[0] + (ur[0]-ll[0])*i/mg[0];
|
|
|
441 n = wdf(x,l0,0);
|
|
|
442
|
|
|
443 t = sqrt(innerprod(l0,l0,n));
|
|
|
444 for (j=0; j<n; j++) l0[j] /= t;
|
|
|
445
|
|
|
446 if (i>0)
|
|
|
447 { t = 0.0;
|
|
|
448 for (j=0; j<n; j++) t += (l1[j]-l0[j])*(l1[j]-l0[j]);
|
|
|
449 sum += 2*asin(sqrt(t)/2);
|
|
|
450 }
|
|
|
451 }
|
|
|
452 kap[0] = sum;
|
|
|
453 if (kap_terms<=1) return(1);
|
|
|
454 kap[1] = 0.0;
|
|
|
455 lap[0] = 2.0;
|
|
|
456 return(2);
|
|
|
457 }
|
|
|
458
|
|
|
459 int tube_constants(f,d,m,ev,mg,fl,kap,wk,terms,uc)
|
|
|
460 double *fl, *kap, *wk;
|
|
|
461 int d, m, ev, *mg, (*f)(), terms, uc;
|
|
|
462 { int aw, deb=0;
|
|
|
463 double k0[4], l0[3], m0[2], n0[1], z[TUBE_MXDIM];
|
|
|
464
|
|
|
465 wdf = f;
|
|
|
466
|
|
|
467 aw = (wk==NULL);
|
|
|
468 if (aw) {
|
|
|
469 wk = (double *)calloc(k0_reqd(d,m,uc),sizeof(double));
|
|
|
470 if ( wk == NULL ) {
|
|
|
471 printf("Problem allocating memory for wk\n");fflush(stdout);
|
|
|
472 }
|
|
|
473 }
|
|
|
474 assignk0(wk,d,m);
|
|
|
475
|
|
|
476 k0[0] = k0[1] = k0[2] = k0[3] = 0.0;
|
|
|
477 l0[0] = l0[1] = l0[2] = 0.0;
|
|
|
478 m0[0] = m0[1] = 0.0;
|
|
|
479 n0[0] = 0.0;
|
|
|
480
|
|
|
481 use_covar = uc;
|
|
|
482 kap_terms = terms;
|
|
|
483 if ((kap_terms <=0) | (kap_terms >= 5))
|
|
|
484 mut_printf("Warning: terms = %2d\n",kap_terms);
|
|
|
485
|
|
|
486 switch(ev)
|
|
|
487 {
|
|
|
488 case IMONTE:
|
|
|
489 monte(k0x,fl,&fl[d],d,k0,mg[0]);
|
|
|
490 break;
|
|
|
491 case ISPHERIC:
|
|
|
492 if (d==2) integ_disc(k0x,l1x,fl,k0,l0,mg);
|
|
|
493 if (d==3) integ_sphere(k0x,l1x,fl,k0,l0,mg);
|
|
|
494 break;
|
|
|
495 case ISIMPSON:
|
|
|
496 if (use_covar) simpson4(k0x,l1x,m0x,n0x,fl,&fl[d],d,k0,l0,m0,n0,mg,z);
|
|
|
497 else simpson4(k0x,l1x,m0x,n0x,fl,&fl[d],d,k0,l0,m0,n0,mg,z);
|
|
|
498 break;
|
|
|
499 case IDERFREE:
|
|
|
500 kodf(fl,&fl[d],mg,k0,l0);
|
|
|
501 break;
|
|
|
502 default:
|
|
|
503 mut_printf("Unknown integration type in tube_constants().\n");
|
|
|
504 }
|
|
|
505
|
|
|
506 if (deb>0)
|
|
|
507 { mut_printf("constants:\n");
|
|
|
508 mut_printf(" k0: %8.5f %8.5f %8.5f %8.5f\n",k0[0],k0[1],k0[2],k0[3]);
|
|
|
509 mut_printf(" l0: %8.5f %8.5f %8.5f\n",l0[0],l0[1],l0[2]);
|
|
|
510 mut_printf(" m0: %8.5f %8.5f\n",m0[0],m0[1]);
|
|
|
511 mut_printf(" n0: %8.5f\n",n0[0]);
|
|
|
512 if (d==2) mut_printf(" check: %8.5f\n",(k0[0]+k0[2]+l0[1]+m0[0])/(2*PI));
|
|
|
513 if (d==3) mut_printf(" check: %8.5f\n",(l0[0]+l0[2]+m0[1]+n0[0])/(4*PI));
|
|
|
514 }
|
|
|
515
|
|
|
516 if (aw) free(wk);
|
|
|
517
|
|
|
518 kap[0] = k0[0];
|
|
|
519 if (kap_terms==1) return(1);
|
|
|
520 kap[1] = l0[0]/2;
|
|
|
521 if ((kap_terms==2) | (d==1)) return(2);
|
|
|
522 kap[2] = (k0[2]+l0[1]+m0[0])/(2*PI);
|
|
|
523 if ((kap_terms==3) | (d==2)) return(3);
|
|
|
524 kap[3] = (l0[2]+m0[1]+n0[0])/(4*PI);
|
|
|
525 return(4);
|
|
|
526 }
|
|
|
527 /*
|
|
|
528 * Copyright 1996-2006 Catherine Loader.
|
|
|
529 */
|
|
|
530 /*
|
|
|
531 * Copyright (c) 1998-2006 Catherine Loader.
|
|
|
532 *
|
|
|
533 * Computes the critical values from constants kappa0 etc
|
|
|
534 * and significance level.
|
|
|
535 */
|
|
|
536
|
|
|
537 #include <math.h>
|
|
|
538 #include "tube.h"
|
|
|
539
|
|
|
540 #define LOGPI 1.144729885849400174143427
|
|
|
541
|
|
|
542 /* area(d) = 2 pi^(d/2) / Gamma(d/2)
|
|
|
543 * = surface area of unit sphere in R^d
|
|
|
544 */
|
|
|
545 static double A[10] =
|
|
|
546 { 1, /* d=0, whatever */
|
|
|
547 2,
|
|
|
548 6.2831853071795864770, /* 2*pi */
|
|
|
549 12.566370614359172954, /* 4*pi */
|
|
|
550 19.739208802178717238, /* 2*pi^2 */
|
|
|
551 26.318945069571622985, /* 8/3*pi^2 */
|
|
|
552 31.006276680299820177, /* pi^3 */
|
|
|
553 33.073361792319808190, /* 16/15*pi^3 */
|
|
|
554 32.469697011334145747, /* 1/3*pi^4 */
|
|
|
555 29.686580124648361825 /* 32/105*pi^4 */
|
|
|
556 };
|
|
|
557
|
|
|
558 double area(d)
|
|
|
559 int d;
|
|
|
560 { if (d<10) return(A[d]);
|
|
|
561 return(2*exp(d*LOGPI/2.0-mut_lgammai(d)));
|
|
|
562 }
|
|
|
563
|
|
|
564 double tailp_uniform(c,k0,m,d,s,n)
|
|
|
565 double c, *k0, n;
|
|
|
566 int m, d, s;
|
|
|
567 { int i;
|
|
|
568 double p;
|
|
|
569 p = 0.0;
|
|
|
570 for (i=0; i<m; i++) if (k0[i] != 0.0)
|
|
|
571 p += k0[i] * ibeta(1-c*c,(n-d+i-1)/2.0,(d+1-i)/2.0) / area(d+1-i);
|
|
|
572 return( (s==TWO_SIDED) ? 2*p : p );
|
|
|
573 }
|
|
|
574
|
|
|
575 double tailp_gaussian(c,k0,m,d,s,n)
|
|
|
576 double c, *k0, n;
|
|
|
577 int m, d, s;
|
|
|
578 { int i;
|
|
|
579 double p;
|
|
|
580 p = 0.0;
|
|
|
581 for (i=0; i<m; i++) if (k0[i] != 0.0)
|
|
|
582 p += k0[i] * (1-pchisq(c*c,(double) d+1-i)) / area(d+1-i);
|
|
|
583 return( (s==TWO_SIDED) ? 2*p : p );
|
|
|
584 }
|
|
|
585
|
|
|
586 double tailp_tprocess(c,k0,m,d,s,n)
|
|
|
587 double c, *k0, n;
|
|
|
588 int m, d, s;
|
|
|
589 { int i;
|
|
|
590 double p;
|
|
|
591 p = 0.0;
|
|
|
592 for (i=0; i<m; i++) if (k0[i] != 0.0)
|
|
|
593 p += k0[i] * (1-pf(c*c/(d+1-i),(double) d+1-i, n)) / area(d+1-i);
|
|
|
594 return( (s==TWO_SIDED) ? 2*p : p );
|
|
|
595 }
|
|
|
596
|
|
|
597 double taild_uniform(c,k0,m,d,s,n)
|
|
|
598 double c, *k0, n;
|
|
|
599 int m, d, s;
|
|
|
600 { int i;
|
|
|
601 double p;
|
|
|
602 p = 0.0;
|
|
|
603 for (i=0; i<m; i++) if (k0[i] != 0.0)
|
|
|
604 p += k0[i] * 2*c*dbeta(1-c*c,(n-d+i-1)/2.0,(d+1-i)/2.0,0) / area(d+1-i);
|
|
|
605 return( (s==TWO_SIDED) ? 2*p : p );
|
|
|
606 }
|
|
|
607
|
|
|
608 double taild_gaussian(c,k0,m,d,s,n)
|
|
|
609 double c, *k0, n;
|
|
|
610 int m, d, s;
|
|
|
611 { int i;
|
|
|
612 double p;
|
|
|
613 p = 0.0;
|
|
|
614 for (i=0; i<m; i++) if (k0[i] != 0.0)
|
|
|
615 p += k0[i] * 2*c*dchisq(c*c,(double) d+1-i,0) / area(d+1-i);
|
|
|
616 return( (s==TWO_SIDED) ? 2*p : p );
|
|
|
617 }
|
|
|
618
|
|
|
619 double taild_tprocess(c,k0,m,d,s,n)
|
|
|
620 double c, *k0, n;
|
|
|
621 int m, d, s;
|
|
|
622 { int i;
|
|
|
623 double p;
|
|
|
624 p = 0.0;
|
|
|
625 for (i=0; i<m; i++) if (k0[i] != 0.0)
|
|
|
626 p += k0[i] * 2*c*df(c*c/(d+1-i),(double) d+1-i, n,0) / ((d+1-i)*area(d+1-i));
|
|
|
627 return( (s==TWO_SIDED) ? 2*p : p );
|
|
|
628 }
|
|
|
629
|
|
|
630 double tailp(c,k0,m,d,s,nu, process)
|
|
|
631 double c, *k0, nu;
|
|
|
632 int m, d, s, process;
|
|
|
633 { switch(process)
|
|
|
634 { case UNIF: return(tailp_uniform(c,k0,m,d,s,nu));
|
|
|
635 case GAUSS: return(tailp_gaussian(c,k0,m,d,s,nu));
|
|
|
636 case TPROC: return(tailp_tprocess(c,k0,m,d,s,nu));
|
|
|
637 }
|
|
|
638 mut_printf("taild: unknown process.\n");
|
|
|
639 return(0.0);
|
|
|
640 }
|
|
|
641
|
|
|
642 double taild(c,k0,m,d,s,nu, process)
|
|
|
643 double c, *k0, nu;
|
|
|
644 int m, d, s, process;
|
|
|
645 { switch(process)
|
|
|
646 { case UNIF: return(taild_uniform(c,k0,m,d,s,nu));
|
|
|
647 case GAUSS: return(taild_gaussian(c,k0,m,d,s,nu));
|
|
|
648 case TPROC: return(taild_tprocess(c,k0,m,d,s,nu));
|
|
|
649 }
|
|
|
650 mut_printf("taild: unknown process.\n");
|
|
|
651 return(0.0);
|
|
|
652 }
|
|
|
653
|
|
|
654 double critval(alpha,k0,m,d,s,nu,process)
|
|
|
655 double alpha, *k0, nu;
|
|
|
656 int m, d, s, process;
|
|
|
657 { double c, cn, c0, c1, tp, td;
|
|
|
658 int j, maxit;
|
|
|
659 double (*tpf)(), (*tdf)();
|
|
|
660
|
|
|
661 maxit = 20;
|
|
|
662 if (m<0)
|
|
|
663 { mut_printf("critval: no terms?\n");
|
|
|
664 return(2.0);
|
|
|
665 }
|
|
|
666 if (m>d+1) m = d+1;
|
|
|
667 if ((alpha<=0) | (alpha>=1))
|
|
|
668 { mut_printf("critval: invalid alpha %8.5f\n",alpha);
|
|
|
669 return(2.0);
|
|
|
670 }
|
|
|
671 if (alpha>0.5)
|
|
|
672 mut_printf("critval: A mighty large tail probability alpha=%8.5f\n",alpha);
|
|
|
673 if (m==0) { d = 0; k0[0] = 1; m = 1; }
|
|
|
674
|
|
|
675 switch(process)
|
|
|
676 { case UNIF:
|
|
|
677 c = 0.5; c0 = 0.0; c1 = 1.0;
|
|
|
678 tpf = tailp_uniform;
|
|
|
679 tdf = taild_uniform;
|
|
|
680 break;
|
|
|
681 case GAUSS:
|
|
|
682 c = 2.0; c0 = 0.0; c1 = 0.0;
|
|
|
683 tpf = tailp_gaussian;
|
|
|
684 tdf = taild_gaussian;
|
|
|
685 break;
|
|
|
686 case TPROC:
|
|
|
687 c = 2.0; c0 = 0.0; c1 = 0.0;
|
|
|
688 tpf = tailp_tprocess;
|
|
|
689 tdf = taild_tprocess;
|
|
|
690 break;
|
|
|
691 default:
|
|
|
692 mut_printf("critval: unknown process.\n");
|
|
|
693 return(0.0);
|
|
|
694 }
|
|
|
695
|
|
|
696 for (j=0; j<maxit; j++)
|
|
|
697 { tp = tpf(c,k0,m,d,s,nu)-alpha;
|
|
|
698 td = tdf(c,k0,m,d,s,nu);
|
|
|
699 if (tp>0) c0 = c;
|
|
|
700 if (tp<0) c1 = c;
|
|
|
701 cn = c + tp/td;
|
|
|
702 if (cn<c0) cn = (c+c0)/2;
|
|
|
703 if ((c1>0.0) && (cn>c1)) cn = (c+c1)/2;
|
|
|
704 c = cn;
|
|
|
705 if (fabs(tp/alpha)<1.0e-10) return(c);
|
|
|
706 }
|
|
|
707 return(c);
|
|
|
708 }
|
