annotate rDiff/src/locfit/Source/libtube.c @ 3:29a698dc5c7e default tip

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