SCM

SCM Repository

[matrix] Annotation of /pkg/src/ssclme.c
ViewVC logotype

Annotation of /pkg/src/ssclme.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 21 - (view) (download) (as text)

1 : bates 10 #include "ssclme.h"
2 :    
3 :     /**
4 :     * Check for a nested series of grouping factors in the sparse,
5 :     * symmetric representation of the pairwise cross-tabulations.
6 :     *
7 :     * @param n size of pairwise cross-tabulation matrix
8 :     * @param nf number of groups of columns in pairwise cross-tabulation
9 :     * @param upper non-zero if the upper triangle is stored
10 :     * @param Ap array of pointers to columns
11 :     * @param Ai row indices
12 :     * @param Gp array of pointers to groups
13 :     *
14 :     * @return 0 for non-nested groups, 1 for nested groups
15 :     */
16 :     static
17 :     int ctab_isNested(int n, int nf, int upper,
18 :     const int Ap[], const int Ai[], const int Gp[])
19 :     {
20 :     if (nf > 1) { /* single factor always nested */
21 :     int i;
22 :     if (upper) {
23 :     int *nnz = (int *) R_alloc(n, sizeof(int)), nz = Ap[n];
24 :     /* count number of nonzeros in each row */
25 :     for (i = 0; i < n; i++) nnz[i] = 0;
26 :     for (i = 0; i < nz; i++) nnz[Ai[i]]++;
27 :     for (i = 0; i < nf; i++) {
28 :     int j, p2 = Gp[i+1], target = nf - i;
29 :     for (j = Gp[i]; j < p2; j++) {
30 :     if (nnz[j] != target) return 0;
31 :     }
32 :     }
33 :     } else { /* lower triangle - the easy case */
34 :     for (i = 0; i < nf; i++) {
35 :     int j, p2 = Gp[i+1], target = nf - i;
36 :     for (j = Gp[i]; j < p2; j++) {
37 :     if ((Ap[j+1] - Ap[j]) != target)
38 :     return 0;
39 :     }
40 :     }
41 :     }
42 :     }
43 :     return 1;
44 :     }
45 :    
46 :     /**
47 :     * Determine if a fill-reducing permutation is needed for the pairwise
48 :     * cross-tabulation matrix. If so, determine such a permutation
49 :     * (using Metis) then separate the groups.
50 :     *
51 :     * @param ctab pointer to a pairwise cross-tabulation object
52 :     *
53 :     * @return pointer to an integer R vector.
54 :     */
55 :     static
56 :     SEXP ctab_permute(SEXP ctab)
57 :     {
58 :     SEXP val, GpSl = GET_SLOT(ctab, Matrix_GpSym);
59 :     int *Ai = INTEGER(GET_SLOT(ctab, Matrix_iSym)),
60 :     *Ap = INTEGER(GET_SLOT(ctab, Matrix_pSym)),
61 :     *Gp = INTEGER(GpSl),
62 :     *perm,
63 :     *work,
64 :     i,
65 :     j,
66 :     n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
67 :     nf = length(GpSl) - 1,
68 :     nz = Ap[n], /* number of non-zeros */
69 :     pos;
70 :    
71 :     if (ctab_isNested(n, nf, 1, Ap, Ai, Gp))
72 :     return allocVector(INTSXP, 0);
73 :     val = allocVector(INTSXP, n);
74 :     perm = INTEGER(val);
75 :     work = (int *) R_alloc(n, sizeof(int));
76 :     ssc_metis_order(n, nz, Ap, Ai, work, perm); /* perm gets inverse perm */
77 :     /* work now contains desired permutation but with groups scrambled */
78 :    
79 :     /* copy work into perm preserving the order of the groups */
80 :     pos = 0; /* position in new permutation */
81 :     for (i = 0; i < nf; i++) {
82 :     for (j = 0; j < n; j++) {
83 :     int jj = work[j];
84 :     if (Gp[i] <= jj && jj < Gp[i+1]) {
85 :     perm[pos] = jj;
86 :     pos++;
87 :     }
88 :     }
89 :     }
90 :     return val;
91 :     }
92 :    
93 :     static
94 :     void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)
95 :     {
96 :     int *snc, i, copyonly = 1;
97 :    
98 :     for (i = 0; i < nf; i++) {
99 :     if (nc[i] > 1) copyonly = 0;
100 :     }
101 :     if (copyonly) {
102 :     SET_SLOT(ssc, Matrix_pSym, duplicate(GET_SLOT(ctab, Matrix_pSym)));
103 :     SET_SLOT(ssc, Matrix_iSym, duplicate(GET_SLOT(ctab, Matrix_iSym)));
104 :     SET_SLOT(ssc, Matrix_xSym, duplicate(GET_SLOT(ctab, Matrix_xSym)));
105 :     SET_SLOT(ssc, Matrix_DimSym,
106 :     duplicate(GET_SLOT(ctab, Matrix_DimSym)));
107 :     SET_SLOT(ssc, Matrix_GpSym, duplicate(GET_SLOT(ctab, Matrix_GpSym)));
108 :     } else {
109 :     int
110 :     *GpIn = INTEGER(GET_SLOT(ctab, Matrix_GpSym)),
111 :     *GpOut,
112 :     *AiIn = INTEGER(GET_SLOT(ctab, Matrix_iSym)),
113 :     *AiOut,
114 :     *ApIn = INTEGER(GET_SLOT(ctab, Matrix_pSym)),
115 :     *ApOut,
116 :     nIn = GpIn[nf], nOut, nzOut,
117 :     *dims,
118 :     *map = Calloc(nIn + 1, int), /* column number map */
119 :     *ncc = Calloc(nIn, int); /* number of columns out for this
120 :     * col in */
121 :    
122 :     SET_SLOT(ssc, Matrix_GpSym, allocVector(INTSXP, nf + 1));
123 :     GpOut = INTEGER(GET_SLOT(ssc, Matrix_GpSym));
124 :     map[0] = GpOut[0] = 0;
125 :     for (i = 0; i < nf; i++) {
126 :     int j, GpIni = GpIn[i], GpInip1 = GpIn[i+1], nci = nc[i];
127 :     GpOut[i+1] = GpOut[i] + (GpInip1 - GpIni) * nci;
128 :     for (j = GpIni; j < GpInip1; j++) {
129 :     ncc[j] = nci;
130 :     map[j+1] = map[j] + nci;
131 :     }
132 :     }
133 :     nOut = GpOut[nf]; /* size of output matrix */
134 :     SET_SLOT(ssc, Matrix_DimSym,
135 :     duplicate(GET_SLOT(ctab, Matrix_DimSym)));
136 :     dims = INTEGER(GET_SLOT(ssc, Matrix_DimSym));
137 :     dims[0] = dims[1] = nOut;
138 :     SET_SLOT(ssc, Matrix_pSym, allocVector(INTSXP, nOut + 1));
139 :     ApOut = INTEGER(GET_SLOT(ssc, Matrix_pSym));
140 :     ApOut[0] = 0;
141 :     for (i = 0; i < nf; i++) { /* determine the column pointers */
142 :     int j, jout = GpOut[i], nci = nc[i], p2 = GpIn[i+1];
143 :     for (j = GpIn[i]; j < p2; j++) {
144 :     int k, nj = 0, p3 = ApIn[j+1];
145 :     for (k = ApIn[j]; k < p3; k++) {
146 :     nj += ncc[AiIn[k]];
147 :     }
148 :     nj -= nci - 1;
149 :     ApOut[jout+1] = ApOut[jout] + nj;
150 :     jout++;
151 :     for (k = 1; k < nci; k++) {
152 :     ApOut[jout+1] = ApOut[jout] + nj + k;
153 :     jout++;
154 :     }
155 :     }
156 :     }
157 :     nzOut = ApOut[nOut]; /* number of non-zeros in output */
158 :     SET_SLOT(ssc, Matrix_xSym, allocVector(REALSXP, nzOut));
159 :     memset(REAL(GET_SLOT(ssc, Matrix_xSym)), 0,
160 :     sizeof(double) * nzOut);
161 :     SET_SLOT(ssc, Matrix_iSym, allocVector(INTSXP, nzOut));
162 :     AiOut = INTEGER(GET_SLOT(ssc, Matrix_iSym));
163 :     for (i = 0; i < nf; i++) { /* fill in the rows */
164 :     int j, jj, nci = nc[i], p2 = GpIn[i+1];
165 :     for (j = GpIn[i]; j < p2; j++) { /* first col for input col */
166 :     int ii = AiIn[j], mj = map[j], ncci = ncc[ii],
167 :     pos = ApOut[mj];
168 :     AiOut[pos++] = map[ii];
169 :     if (ii < j) { /* above the diagonal */
170 :     for (jj = 1; jj < ncci; jj++) {
171 :     AiOut[pos+1] = AiOut[pos] + 1;
172 :     pos++;
173 :     }
174 :     }
175 :     for (jj = 1; jj < nci; jj++) { /* repeat the column adding
176 :     * another diagonal element */
177 :     int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];
178 :     Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);
179 :     AiOut[ApOut[mjj + 1] - 1] = mjj; /* maybe mjj-1? */
180 :     }
181 :     }
182 :     }
183 :     Free(map); Free(ncc);
184 :     }
185 :     SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));
186 :     snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
187 :     for (i = 0; i <= nf; i++) {
188 :     snc[i] = nc[i];
189 :     }
190 :     }
191 :    
192 :     static
193 :     void ssclme_fill_LIp(int n, const int Parent[], int LIp[])
194 :     {
195 :     int *sz = Calloc(n, int), i;
196 :     for (i = n - 1; i >= 0; i--) {
197 :     sz[i] = (Parent[i] < 0) ? 0 : 1 + sz[Parent[i]];
198 :     }
199 :     LIp[0] = 0;
200 :     for (i = 0; i < n; i++) LIp[i+1] = LIp[i] + sz[i];
201 :     Free(sz);
202 :     }
203 :    
204 :     static
205 :     void ssclme_fill_LIi(int n, const int Parent[], const int LIp[], int LIi[])
206 :     {
207 :     int i;
208 :     for (i = n; i > 0; i--) {
209 :     int im1 = i - 1, Par = Parent[im1];
210 :     if (Par >= 0) {
211 :     LIi[LIp[im1]] = Par;
212 :     Memcpy(LIi + LIp[im1] + 1, LIi + LIp[Par],
213 :     LIp[Par + 1] - LIp[Par]);
214 :     }
215 :     }
216 :     }
217 :    
218 :     SEXP
219 :     ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)
220 :     {
221 :     SEXP ctab, nms, ssc, tmp,
222 : bates 21 val = PROTECT(allocVector(VECSXP, 2)),
223 :     dd = PROTECT(allocVector(INTSXP, 3)); /* dimensions of 3-D arrays */
224 : bates 10 int *Ai, *Ap, *Gp, *LIp, *Lp, *Parent,
225 : bates 21 *nc, Lnz, i, nf = length(facs), nzcol, pp1,
226 :     *dims = INTEGER(dd);
227 : bates 10
228 :     if (length(ncv) != (nf + 1))
229 :     error("length of nc (%d) should be length of facs (%d) + 1",
230 :     length(ncv), nf);
231 :     SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS("ssclme")));
232 :     ssc = VECTOR_ELT(val, 0);
233 :     /* Pairwise cross-tabulation */
234 :     ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));
235 :     SET_VECTOR_ELT(val, 1, ctab_permute(ctab));
236 :     if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */
237 :     ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
238 :     1, INTEGER(VECTOR_ELT(val, 1)),
239 :     INTEGER(GET_SLOT(ctab, Matrix_pSym)),
240 :     INTEGER(GET_SLOT(ctab, Matrix_iSym)));
241 :     }
242 :     ssclme_copy_ctab(nf, INTEGER(ncv), ctab, ssc);
243 :     UNPROTECT(1); /* ctab */
244 :    
245 :     nzcol = INTEGER(GET_SLOT(ssc, Matrix_DimSym))[1];
246 :     Gp = INTEGER(GET_SLOT(ssc, Matrix_GpSym));
247 :     Ap = INTEGER(GET_SLOT(ssc, Matrix_pSym));
248 :     Ai = INTEGER(GET_SLOT(ssc, Matrix_iSym));
249 :     nc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
250 :     nc[nf + 1] = length(VECTOR_ELT(facs, 0)); /* number of observations */
251 :     /* Create slots */
252 :     pp1 = nc[nf];
253 :     SET_SLOT(ssc, Matrix_XtXSym, allocMatrix(REALSXP, pp1, pp1));
254 :     SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));
255 :     SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));
256 :     SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));
257 :     /* Zero the symmetric matrices (for cosmetic reasons only). */
258 :     memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,
259 :     sizeof(double) * pp1 * pp1);
260 :     memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,
261 :     sizeof(double) * pp1 * pp1);
262 :     SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));
263 :     Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));
264 :     SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol));
265 :     Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));
266 :     SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));
267 :     SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));
268 :     ldl_symbolic(nzcol, Ap, Ai, Lp, Parent,
269 :     (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */
270 :     (int *) R_alloc(nzcol, sizeof(int)), /* Flag */
271 :     (int *) NULL, (int *) NULL); /* P & Pinv */
272 :     Lnz = Lp[nzcol];
273 :     SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));
274 :     SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));
275 :     SET_SLOT(ssc, Matrix_OmegaSym, allocVector(VECSXP, nf));
276 :     tmp = GET_SLOT(ssc, Matrix_OmegaSym);
277 :     setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));
278 :     for (i = 0; i < nf; i++) {
279 :     SET_VECTOR_ELT(tmp, i, allocMatrix(REALSXP, nc[i], nc[i]));
280 :     memset(REAL(VECTOR_ELT(tmp, i)), 0,
281 :     sizeof(double) * nc[i] * nc[i]);
282 :     }
283 :     SET_SLOT(ssc, Matrix_devianceSym, allocVector(REALSXP, 2));
284 :     tmp = GET_SLOT(ssc, Matrix_devianceSym);
285 :     setAttrib(tmp, R_NamesSymbol, allocVector(STRSXP, 2));
286 :     nms = getAttrib(tmp, R_NamesSymbol);
287 :     SET_STRING_ELT(nms, 0, mkChar("ML"));
288 :     SET_STRING_ELT(nms, 1, mkChar("REML"));
289 :     SET_SLOT(ssc, Matrix_devCompSym, allocVector(REALSXP, 4));
290 :     SET_SLOT(ssc, Matrix_statusSym, allocVector(LGLSXP, 2));
291 :     tmp = GET_SLOT(ssc, Matrix_statusSym);
292 :     LOGICAL(tmp)[0] = LOGICAL(tmp)[1] = 0;
293 :     setAttrib(tmp, R_NamesSymbol, allocVector(STRSXP, 2));
294 :     nms = getAttrib(tmp, R_NamesSymbol);
295 :     SET_STRING_ELT(nms, 0, mkChar("factored"));
296 :     SET_STRING_ELT(nms, 1, mkChar("inverted"));
297 :     SET_SLOT(ssc, Matrix_bVarSym, allocVector(VECSXP, nf));
298 :     tmp = GET_SLOT(ssc, Matrix_bVarSym);
299 :     setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));
300 :     for (i = 0; i < nf; i++) {
301 : bates 21 int nci = nc[i], mi = (Gp[i+1] - Gp[i])/nc[i];
302 : bates 10
303 : bates 21 dims[0] = dims[1] = nci;
304 :     dims[2] = mi;
305 :     SET_VECTOR_ELT(tmp, i, allocArray(REALSXP, dd));
306 : bates 10 memset(REAL(VECTOR_ELT(tmp, i)), 0,
307 : bates 21 sizeof(double) * nci * nci * mi);
308 : bates 10 }
309 :     SET_SLOT(ssc, Matrix_LIpSym, allocVector(INTSXP, nzcol + 1));
310 :     LIp = INTEGER(GET_SLOT(ssc, Matrix_LIpSym));
311 :     ssclme_fill_LIp(nzcol, Parent, LIp);
312 :     if (asInteger(threshold) > (Lnz = LIp[nzcol])) {
313 :     SET_SLOT(ssc, Matrix_LIiSym, allocVector(INTSXP, Lnz));
314 :     ssclme_fill_LIi(nzcol, Parent, LIp,
315 :     INTEGER(GET_SLOT(ssc, Matrix_LIiSym)));
316 :     SET_SLOT(ssc, Matrix_LIxSym, allocVector(REALSXP, Lnz));
317 :     memset(REAL(GET_SLOT(ssc, Matrix_LIxSym)), 0,
318 :     sizeof(double) * Lnz);
319 :     }
320 : bates 21 UNPROTECT(2);
321 : bates 10 return val;
322 :     }
323 :    
324 :     static
325 :     void bVj_to_A(int ncj, int Gpj, int Gpjp, const double bVj[],
326 :     const int Ap[], const int Ai[], double Ax[])
327 :     {
328 :     int i, diag, k;
329 :     for (i = Gpj; i < Gpjp; i += ncj) {
330 :     for (k = 0; k < ncj; k++) {
331 :     diag = Ap[i + k + 1] - 1;
332 :     if (Ai[diag] != i+k)
333 :     error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",
334 :     diag, i+k, Ai[diag]);
335 :     Memcpy(Ax + diag - k, bVj + (i+k-Gpj)*ncj, k + 1);
336 :     }
337 :     }
338 :     }
339 :    
340 :     SEXP
341 :     ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)
342 :     {
343 :     SEXP
344 :     bVar = GET_SLOT(x, Matrix_bVarSym);
345 :     int
346 :     *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
347 :     *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
348 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
349 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
350 :     i, j, k,
351 :     ione = 1,
352 :     nf = length(mmats) - 1,
353 :     nobs = nc[nf + 1],
354 :     nzcol = Gp[nf],
355 :     nz = Ap[nzcol],
356 :     pp1 = nc[nf];
357 :     double
358 :     **Z = Calloc(nf + 1, double *),
359 :     *Ax = REAL(GET_SLOT(x, Matrix_xSym)),
360 :     *XtX = REAL(GET_SLOT(x, Matrix_XtXSym)),
361 :     *ZtX = REAL(GET_SLOT(x, Matrix_ZtXSym)),
362 :     one = 1.0,
363 :     zero = 0.0;
364 :    
365 :     for (i = 0; i <= nf; i++) {
366 :     int *dims = INTEGER(getAttrib(VECTOR_ELT(mmats, i), R_DimSymbol)),
367 :     nci = nc[i];
368 :     if (nobs != dims[0])
369 :     error("Expected %d rows in the %d'th model matrix. Got %d",
370 :     nobs, i+1, dims[0]);
371 :     if (nci != dims[1])
372 :     error("Expected %d columns in the %d'th model matrix. Got %d",
373 :     nci, i+1, dims[1]);
374 :     Z[i] = REAL(VECTOR_ELT(mmats, i));
375 :     }
376 :     /* Create XtX - X is Z[nf] */
377 :     F77_CALL(dsyrk)("U", "T", nc+nf, &nobs, &one,
378 :     Z[nf], &nobs, &zero, XtX, nc + nf);
379 :     /* Zero the accumulators */
380 :     memset((void *) ZtX, 0, sizeof(double) * nzcol * pp1);
381 :     memset((void *) Ax, 0, sizeof(double) * nz);
382 :     for (j = 0; j < nf; j++) { /* Create ZtX */
383 :     int *fpj = INTEGER(VECTOR_ELT(facs, j)), ncj = nc[j],
384 :     Ncj = ncj > 1;
385 :     double
386 :     *bVj = REAL(VECTOR_ELT(bVar, j)),
387 :     *Zj = Z[j],
388 :     *zxj = ZtX + Gp[j];
389 :    
390 :     if (Ncj) { /* bVj will accumulate Z'Z blocks */
391 :     memset(bVj, 0, sizeof(double) * ncj * (Gp[j+1]-Gp[j]));
392 :     }
393 :     for (i = 0; i < nobs; i++) { /* accumulate diagonal of ZtZ */
394 :     int fpji = fpj[i] - 1, /* factor indices are 1-based */
395 :     dind = Ap[Gp[j] + fpji * ncj + 1] - 1;
396 :     if (Ai[dind] != (Gp[j] + fpji * ncj))
397 :     error("logic error in ssclme_update_mm");
398 :     if (Ncj) { /* use bVar to accumulate */
399 :     F77_CALL(dsyrk)("U", "T", &ncj, &ione, &one, Zj+i,
400 :     &nobs, &one, bVj + fpji*ncj*ncj, &ncj);
401 :     } else { /* update scalars directly */
402 :     Ax[dind] += Zj[i] * Zj[i];
403 :     }
404 :     /* update rows of Z'X */
405 :     F77_CALL(dgemm)("T", "N", &ncj, &pp1, &ione, &one,
406 :     Zj + i, &nobs, Z[nf] + i, &nobs,
407 :     &one, zxj + fpji * ncj, &nzcol);
408 :     }
409 :     if (Ncj) bVj_to_A(ncj, Gp[j], Gp[j+1], bVj, Ap, Ai, Ax);
410 :     for (k = j+1; k < nf; k++) { /* off-diagonals */
411 :     int *fpk = INTEGER(VECTOR_ELT(facs, k)),
412 :     *Apk = Ap + Gp[k],
413 :     nck = nc[k];
414 :     double
415 :     *Zk = Z[k];
416 :    
417 :     for (i = 0; i < nobs; i++) {
418 :     int ii, ind = -1, fpji = fpj[i] - 1,
419 :     row = Gp[j] + fpji * ncj,
420 :     fpki = fpk[i] - 1,
421 :     lastind = Apk[fpki + 1];
422 :     for (ii = Apk[fpki]; ii < lastind; ii++) {
423 :     if (Ai[ii] == row) {
424 :     ind = ii;
425 :     break;
426 :     }
427 :     }
428 :     if (ind < 0) error("logic error in ssclme_update_mm");
429 :     if (Ncj || nck > 1) {
430 :     /* FIXME: run a loop to update */
431 :     } else { /* update scalars directly */
432 :     Ax[ind] += Zj[fpji] * Zk[fpki];
433 :     }
434 :     }
435 :     }
436 :     }
437 :     Free(Z);
438 :     return R_NilValue;
439 :     }
440 :    
441 :     SEXP ssclme_inflate_and_factor(SEXP lme)
442 :     {
443 :     SEXP
444 :     GpSlot = GET_SLOT(lme, Matrix_GpSym),
445 :     Omega = GET_SLOT(lme, Matrix_OmegaSym);
446 :     int n = INTEGER(GET_SLOT(lme, Matrix_DimSym))[1];
447 :     int
448 :     *Ai = INTEGER(GET_SLOT(lme, Matrix_iSym)),
449 :     *Ap = INTEGER(GET_SLOT(lme, Matrix_pSym)),
450 :     *Flag = Calloc(n, int),
451 :     *Gp = INTEGER(GpSlot),
452 :     *Lnz = Calloc(n, int),
453 :     *Pattern = Calloc(n, int),
454 :     *nc = INTEGER(GET_SLOT(lme, Matrix_ncSym)),
455 :     j,
456 :     nf = length(GpSlot) - 1;
457 :     double
458 :     *D = REAL(GET_SLOT(lme, Matrix_DSym)),
459 :     *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),
460 :     *Y = Calloc(n, double),
461 :     *xcp = Calloc(Ap[n], double);
462 :    
463 :     Memcpy(xcp, REAL(GET_SLOT(lme, Matrix_xSym)), Ap[n]);
464 :     for (j = 0; j < nf; j++) {
465 :     int diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];
466 :     double *omgj = REAL(VECTOR_ELT(Omega, j));
467 :    
468 :     for (i = Gp[j]; i < G2; i += ncj) {
469 :     for (k = 0; k < ncj; k++) {
470 :     diag = Ap[i + k + 1] - 1;
471 :     if (Ai[diag] != i+k)
472 :     error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",
473 :     diag, i+k, Ai[diag]);
474 :     for (ii = 0; ii <= k; ii++) {
475 :     xcp[diag + ii - k] += omgj[k*ncj + ii];
476 :     }
477 :     }
478 :     }
479 :     }
480 :     j = ldl_numeric(n, Ap, Ai, xcp,
481 :     INTEGER(GET_SLOT(lme, Matrix_LpSym)),
482 :     INTEGER(GET_SLOT(lme, Matrix_ParentSym)),
483 :     Lnz, INTEGER(GET_SLOT(lme, Matrix_LiSym)),
484 :     REAL(GET_SLOT(lme, Matrix_LxSym)),
485 :     D, Y, Pattern, Flag,
486 :     (int *) NULL, (int *) NULL); /* P & Pinv */
487 :     if (j != n)
488 :     error("rank deficiency of ZtZ+W detected at column %d",
489 :     j + 1);
490 :     for (j = 0; j < n; j++) DIsqrt[j] = 1./sqrt(D[j]);
491 :     Free(Lnz); Free(Flag); Free(Pattern); Free(Y); Free(xcp);
492 :     return R_NilValue;
493 :     }
494 :    
495 :     SEXP ssclme_factor(SEXP lme)
496 :     {
497 :     int *status = LOGICAL(GET_SLOT(lme, Matrix_statusSym));
498 :    
499 :     if (!status[0]) {
500 :     SEXP
501 :     GpSlot = GET_SLOT(lme, Matrix_GpSym),
502 :     Omega = GET_SLOT(lme, Matrix_OmegaSym);
503 :     int
504 :     *Gp = INTEGER(GpSlot),
505 :     *Li = INTEGER(GET_SLOT(lme, Matrix_LiSym)),
506 :     *Lp = INTEGER(GET_SLOT(lme, Matrix_LpSym)),
507 :     *nc = INTEGER(GET_SLOT(lme, Matrix_ncSym)),
508 :     i,
509 :     n = INTEGER(GET_SLOT(lme, Matrix_DimSym))[1],
510 :     nf = length(GpSlot) - 1,
511 :     nobs = nc[nf + 1],
512 :     nreml = nobs + 1 - nc[nf],
513 :     pp1 = nc[nf],
514 :     pp2 = pp1 + 1;
515 :     double
516 :     *D = REAL(GET_SLOT(lme, Matrix_DSym)),
517 :     *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),
518 :     *Lx = REAL(GET_SLOT(lme, Matrix_LxSym)),
519 :     *RXX = REAL(GET_SLOT(lme, Matrix_RXXSym)),
520 :     *RZX = REAL(GET_SLOT(lme, Matrix_RZXSym)),
521 :     *dcmp = REAL(getAttrib(lme, Matrix_devCompSym)),
522 :     *deviance = REAL(getAttrib(lme, Matrix_devianceSym)),
523 :     minus1 = -1.,
524 :     one = 1.;
525 :    
526 :     ssclme_inflate_and_factor(lme);
527 :     /* Accumulate logdet of ZtZ+W */
528 :     dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;
529 :     for (i = 0; i < n; i++) dcmp[0] += log(D[i]);
530 :     /* Accumulate logdet of W */
531 :     for (i = 0; i < nf; i++) {
532 :     int nci = nc[i],
533 :     mi = (Gp[i+1] - Gp[i])/nci;
534 :    
535 :     if (nci < 2) {
536 :     dcmp[1] += mi * log(REAL(VECTOR_ELT(Omega, i))[0]);
537 :     } else {
538 :     int j;
539 :     double
540 :     *tmp = Calloc(nci * nci, double),
541 :     accum = 0.;
542 :     F77_CALL(dpotrf)("U", &nci,
543 :     Memcpy(tmp, REAL(VECTOR_ELT(Omega, i)),
544 :     nci * nci),
545 :     &nci, &j);
546 :     if (j)
547 :     error("Omega[%d] is not positive definite", i + 1);
548 :     for (j = 0; j < nci; j++) {
549 :     accum += 2 * log(tmp[j * (nci + 1)]);
550 :     }
551 :     dcmp[1] += mi * accum;
552 :     Free(tmp);
553 :     }
554 :     }
555 :     /* ldl_lsolve on Z'X */
556 :     Memcpy(RZX, REAL(GET_SLOT(lme, Matrix_ZtXSym)), n * pp1);
557 :     for (i = 0; i < pp1; i++) {
558 :     int j;
559 :     double *RZXi = RZX + i * n;
560 :     ldl_lsolve(n, RZXi, Lp, Li, Lx);
561 :     for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
562 :     }
563 :     /* downdate and factor X'X */
564 :     Memcpy(RXX, REAL(GET_SLOT(lme, Matrix_XtXSym)), pp1 * pp1);
565 :     F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
566 :     RZX, &n, &one, RXX, &pp1);
567 :     F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
568 :     if (i)
569 :     error("DPOTRF returned error code %d", i);
570 :     /* logdet of RXX */
571 :     for (i = 0; i < (pp1 - 1); i++)
572 :     dcmp[2] += 2 * log(RXX[i*pp2]);
573 :     /* logdet of Ryy */
574 :     dcmp[3] = 2. * log(RXX[pp1 * pp1 - 1]);
575 :     deviance[0] = /* ML criterion */
576 :     dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));
577 :     deviance[1] = dcmp[0] - dcmp[1] + /* REML */
578 :     dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));
579 :     status[0] = 1;
580 :     status[1] = 0;
581 :     }
582 :     return R_NilValue;
583 :     }
584 :    
585 :     static
586 :     int ldl_update_ind(int probe, int start, const int ind[])
587 :     {
588 :     while (ind[start] < probe) start++;
589 :     if (ind[start] > probe) error("logic error in ldl_inverse");
590 :     return start;
591 :     }
592 :    
593 :     /**
594 :     * Create the inverse of L and update the diagonal blocks of the inverse
595 :     * of LDL' (=Z'Z+W)
596 :     *
597 :     * @param x pointer to an ssclme object
598 :     *
599 :     * @return R_NilValue (x is updated in place)
600 :    
601 :     */
602 :     SEXP ldl_inverse(SEXP x)
603 :     {
604 :     SEXP
605 :     Gpsl = GET_SLOT(x, Matrix_GpSym),
606 :     LIisl = GET_SLOT(x, Matrix_LIiSym),
607 :     LIpsl = GET_SLOT(x, Matrix_LIpSym),
608 :     bVar = GET_SLOT(x, Matrix_bVarSym);
609 :     int *Gp = INTEGER(Gpsl),
610 :     *Li,
611 :     *LIp = INTEGER(LIpsl), *Lp,
612 :     i,
613 :     nf = length(Gpsl) - 1,
614 :     nzc = length(LIpsl) - 1;
615 :     double
616 :     *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
617 :     *Lx;
618 :    
619 :     ssclme_factor(x);
620 :     if (LIp[nzc] == 0) { /* L and LI are the identity */
621 :     for (i = 0; i < nf; i++) {
622 :     Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],
623 :     Gp[i+1] - Gp[i]);
624 :     }
625 :     return R_NilValue;
626 :     }
627 :     Lp = INTEGER(GET_SLOT(x, Matrix_LpSym));
628 :     Li = INTEGER(GET_SLOT(x, Matrix_LiSym));
629 :     Lx = REAL(GET_SLOT(x, Matrix_LxSym));
630 :     if (length(LIisl) == LIp[nzc]) { /* LIi is filled */
631 :     int *LIi = INTEGER(LIisl),
632 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
633 :     j, jj, k, kk, p1, p2, pi1, pi2;
634 :    
635 :     double *LIx = REAL(GET_SLOT(x, Matrix_LIxSym)),
636 :     one = 1., zero = 0.;
637 :    
638 :     memset(LIx, 0, sizeof(double) * LIp[nzc]);
639 :     /* calculate inverse */
640 :     for (i = 0; i < nzc; i++) {
641 :     p1 = Lp[i]; p2 = Lp[i+1]; pi1 = LIp[i]; pi2 = LIp[i+1];
642 :     /* initialize from unit diagonal term */
643 :     kk = pi1;
644 :     for (j = p1; j < p2; j++) {
645 :     k = Li[j];
646 :     while (LIi[kk] < k && kk < pi2) kk++;
647 :     if (LIi[kk] != k) error("logic error in ldl_inverse");
648 :     LIx[kk] = -Lx[j];
649 :     }
650 :     for (j = pi1; j < pi2; j++) {
651 :     jj = LIi[j];
652 :     p1 = Lp[jj]; p2 = Lp[jj+1];
653 :     kk = j;
654 :     for (jj = p1; jj < p2; jj++) {
655 :     k = Li[jj];
656 :     while (LIi[kk] < k && kk < pi2) kk++;
657 :     if (LIi[kk] != k) error("logic error in ldl_inverse");
658 :     LIx[kk] -= Lx[jj]*LIx[j];
659 :     }
660 :     }
661 :     }
662 :     for (i = 0; i < nf; i++) { /* accumulate bVar */
663 :     int G1 = Gp[i], G2 = Gp[i+1], j, k, kk,
664 :     nci = nc[i], nr, nr1, rr;
665 :     double *bVi = REAL(VECTOR_ELT(bVar, i)), *tmp;
666 :    
667 :     nr = -1;
668 :     for (j = G1; j < G2; j += nci) {
669 :     rr = 1 + LIp[j + 1] - LIp[j];
670 :     if (rr > nr) nr = rr;
671 :     }
672 :     tmp = Calloc(nr * nci, double); /* scratch storage */
673 :     nr1 = nr + 1;
674 :     /* initialize bVi to zero (cosmetic) */
675 :     memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
676 :     for (j = G1; j < G2; j += nci) {
677 :     memset(tmp, 0, sizeof(double) * nr * nci);
678 :     rr = 1 + LIp[j + 1] - LIp[j];
679 :     for (k = 0; k < nci; k++) { /* copy columns */
680 :     tmp[k * nr1] = 1.; /* (unstored) diagonal elt */
681 :     Memcpy(tmp + k*nr1 + 1, LIx + LIp[j + k], rr - k - 1);
682 :     }
683 :     /* scale the rows */
684 :     tmp[0] = DIsqrt[j]; /* first row only has one non-zero */
685 :     for (kk = 1; kk < rr; kk++) {
686 :     for (k = 0; k < nci; k++) {
687 :     tmp[k * nr + kk] *= DIsqrt[LIi[LIp[j] + kk - 1]];
688 :     }
689 :     }
690 :     F77_CALL(dsyrk)("U", "T", &nci, &rr, &one, tmp, &nr,
691 :     &zero, bVi + (j - G1) * nci, &nci);
692 :     F77_CALL(dpotrf)("U", &nci, bVi + (j - G1) * nci,
693 :     &nci, &kk);
694 :     if (kk) /* should never happen */
695 :     error(
696 :     "Rank deficient variance matrix at group %d, level %d",
697 :     i + 1, j + 1);
698 :     }
699 :     }
700 :     return R_NilValue;
701 :     }
702 :     if (length(LIisl)) error("logic error in ssclme_ldl_inverse");
703 :     else { /* LIi and LIx are too big and not used */
704 :     int *counts = Calloc(nzc, int), info, maxod = -1;
705 :     int *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym));
706 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));
707 :     double one = 1.0, zero = 0.;
708 :     /* determine maximum # of off-diagonals */
709 :     for (i = nzc - 1; i >= 0; i--) { /* in a column of L^{-1} */
710 :     counts[i] = (Parent[i] < 0) ? 0 : 1 + counts[Parent[i]];
711 :     if (counts[i] > maxod) maxod = counts[i];
712 :     }
713 :     Free(counts);
714 :    
715 :     for (i = 0; i < nf; i++) {
716 :     int j, jj, k, kk, nci = nc[i], nr, p, p2, pp,
717 :     m = maxod + nci,
718 :     *ind = Calloc(m, int);
719 :     double
720 :     *tmp = Calloc(m * nci, double),
721 :     *mpt = REAL(VECTOR_ELT(bVar, i));
722 :    
723 :     for (j = Gp[i]; j < Gp[i+1]; j += nci) {
724 :     memset((void *) tmp, 0, sizeof(double) * m * nci);
725 :    
726 :     kk = 0; /* ind holds indices of non-zeros */
727 :     jj = j; /* in this block of columns */
728 :     while (jj >= 0) {
729 :     ind[kk++] = jj;
730 :     jj = Parent[jj];
731 :     }
732 :     nr = kk; /* number of non-zeros in this block */
733 :     while (kk < m) ind[kk++] = nzc; /* placeholders */
734 :    
735 :     for (k = 0; k < nci; k++) {
736 :     double *ccol = tmp + k * nr;
737 :    
738 :     ccol[k] = 1.;
739 :     kk = k;
740 :     for (jj = j + k; jj >= 0; jj = Parent[jj]) {
741 :     p2 = Lp[jj+1];
742 :     pp = kk;
743 :     for (p = Lp[jj]; p < p2; p++) {
744 :     pp = ldl_update_ind(Li[p], pp, ind);
745 :     ccol[pp] -= Lx[p] * ccol[kk];
746 :     }
747 :     }
748 :     }
749 :     /* scale rows */
750 :     for (kk = 0; kk < nr; kk++) {
751 :     for (k = 0; k < nci; k++) {
752 :     tmp[k * nr + kk] *= DIsqrt[ind[kk]];
753 :     }
754 :     }
755 :     F77_CALL(dsyrk)("U", "T", &nci, &nr, &one, tmp, &nr,
756 :     &zero, mpt + (j - Gp[i])*nci, &nci);
757 :     F77_CALL(dpotrf)("U", &nci, mpt + (j - Gp[i])*nci,
758 :     &nci, &info);
759 :     if (info) /* should never happen */
760 :     error(
761 :     "Rank deficient variance matrix at group %d, level %d",
762 :     i + 1, j + 1);
763 :     }
764 :     Free(tmp); Free(ind);
765 :     }
766 :     }
767 :     return R_NilValue;
768 :     }
769 :    
770 :     SEXP ssclme_invert(SEXP lme)
771 :     {
772 :     int *status = LOGICAL(GET_SLOT(lme, Matrix_statusSym));
773 :     if (!status[0]) ssclme_factor(lme);
774 :     if (!status[1]) {
775 :     SEXP
776 :     RZXsl = GET_SLOT(lme, Matrix_RZXSym);
777 :     int
778 :     *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
779 :     *Li = INTEGER(GET_SLOT(lme, Matrix_LiSym)),
780 :     *Lp = INTEGER(GET_SLOT(lme, Matrix_LpSym)),
781 :     i,
782 :     n = dims[0],
783 :     pp1 = dims[1];
784 :     double
785 :     *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),
786 :     *Lx = REAL(GET_SLOT(lme, Matrix_LxSym)),
787 :     *RXX = REAL(GET_SLOT(lme, Matrix_RXXSym)),
788 :     *RZX = REAL(RZXsl),
789 :     one = 1.;
790 :    
791 :     F77_CALL(dtrtri)("U", "N", &pp1, RXX, &pp1, &i);
792 :     if (i)
793 :     error("DTRTRI returned error code %d", i);
794 :     F77_CALL(dtrmm)("R", "U", "N", "N", &n, &pp1, &one,
795 :     RXX, &pp1, RZX, &n);
796 :     for (i = 0; i < pp1; i++) {
797 :     int j; double *RZXi = RZX + i * n;
798 :     for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
799 :     ldl_ltsolve(n, RZXi, Lp, Li, Lx);
800 :     }
801 :     ldl_inverse(lme);
802 :     status[1] = 1;
803 :     }
804 :     return R_NilValue;
805 :     }
806 :    
807 :     SEXP ssclme_initial(SEXP x)
808 :     {
809 :     SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),
810 :     Omg = GET_SLOT(x, Matrix_OmegaSym);
811 :     int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
812 :     *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
813 :     *Gp = INTEGER(Gpsl),
814 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
815 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
816 :     i, nf = length(Gpsl) - 1;
817 :    
818 :     double *Ax = REAL(GET_SLOT(x, Matrix_xSym));
819 :    
820 :     for (i = 0; i < nf; i++) {
821 :     int
822 :     Gpi = Gp[i],
823 :     j, k,
824 :     nci = nc[i],
825 :     ncip1 = nci + 1,
826 :     p2 = Gp[i+1];
827 :     double
828 :     mi = 0.375 * ((double) nci)/((double) (p2 - Gpi)),
829 :     *mm = REAL(VECTOR_ELT(Omg, i));
830 :    
831 :     memset((void *) mm, 0, sizeof(double) * nci * nci);
832 :     for (j = Gpi; j < p2; j += nci) {
833 :     for (k = 0; k < nci; k++) {
834 :     int jk = j+k, jj = Ap[jk+1] - 1;
835 :     if (Ai[jj] != jk) error("malformed ZtZ structure");
836 :     mm[k * ncip1] += Ax[jj] * mi;
837 :     }
838 :     }
839 :     }
840 :     status[0] = status[1] = 0;
841 :     return R_NilValue;
842 :     }
843 :    
844 :     /**
845 :     * Extract the conditional estimates of the fixed effects
846 :     * FIXME: Add names
847 :     *
848 :     * @param x Pointer to an ssclme object
849 :     *
850 :     * @return a numeric vector containing the conditional estimates of
851 :     * the fixed effects
852 :     */
853 :     SEXP ssclme_fixef(SEXP x)
854 :     {
855 :     SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
856 :     int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1];
857 :     int j, p = pp1 - 1;
858 :     SEXP val = PROTECT(allocVector(REALSXP, p));
859 :     double
860 :     *RXX = REAL(RXXsl),
861 :     *beta = REAL(val),
862 :     nryyinv; /* negative ryy-inverse */
863 :    
864 :     ssclme_invert(x);
865 :     Memcpy(beta, RXX + p * pp1, p);
866 :     nryyinv = -RXX[pp1*pp1 - 1];
867 :     for (j = 0; j < p; j++) beta[j] /= nryyinv;
868 :     UNPROTECT(1);
869 :     return val;
870 :     }
871 :    
872 :     /**
873 :     * Extract the conditional modes of the random effects.
874 :     * FIXME: Change the returned value to be a named list of matrices
875 :     * with dimnames.
876 :     *
877 :     * @param x Pointer to an ssclme object
878 :     *
879 :     * @return a vector containing the conditional modes of the random effects
880 :     */
881 :     SEXP ssclme_ranef(SEXP x)
882 :     {
883 : bates 21 SEXP RZXsl = GET_SLOT(x, Matrix_RZXSym),
884 :     GpSl = GET_SLOT(x, Matrix_GpSym);
885 : bates 10 int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
886 : bates 21 *Gp = INTEGER(GpSl),
887 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
888 :     i, j,
889 : bates 10 n = dims[0],
890 : bates 21 nf = length(GpSl) - 1,
891 : bates 10 pp1 = dims[1],
892 :     p = pp1 - 1;
893 : bates 21 SEXP val = PROTECT(allocVector(VECSXP, nf));
894 : bates 10 double
895 : bates 21 *b = REAL(RZXsl) + n * p,
896 : bates 10 ryyinv; /* ryy-inverse */
897 :    
898 :     ssclme_invert(x);
899 :     ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];
900 : bates 21 for (i = 0; i < nf; i++) {
901 :     int nci = nc[i], Mi = (Gp[i+1] - Gp[i]), mi = Mi/nci;
902 :     double *mm;
903 :    
904 :     SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, mi));
905 :     mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);
906 :     b += nci * mi;
907 :     for (j = 0; j < Mi; j++) mm[j] /= ryyinv;
908 :     }
909 : bates 10 UNPROTECT(1);
910 :     return val;
911 :     }
912 :     /**
913 :     * Extract the ML or REML conditional estimate of sigma
914 :     *
915 :     * @param x pointer to an ssclme object
916 :     * @param REML logical scalar - TRUE if REML estimates are requested
917 :     *
918 :     * @return numeric scalar
919 :     */
920 :     SEXP ssclme_sigma(SEXP x, SEXP REML)
921 :     {
922 :     SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
923 :     int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],
924 :     nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))
925 :     [length(GET_SLOT(x, Matrix_GpSym))];
926 :    
927 :     ssclme_invert(x);
928 :     return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *
929 :     sqrt((double)(asLogical(REML) ?
930 :     nobs + 1 - pp1 : nobs))));
931 :     }
932 :    
933 :     static
934 :     int coef_length(int nf, const int nc[])
935 :     {
936 :     int i, ans = 0;
937 :     for (i = 0; i < nf; i++) ans += (nc[i] * (nc[i] + 1))/2;
938 :     return ans;
939 :     }
940 :    
941 :     /**
942 :     * Extract the upper triangles of the Omega matrices.
943 :     * (These are not in any sense "coefficients" but the extractor is
944 :     * called coef for historical reasons.)
945 :     *
946 :     * @param x pointer to an ssclme object
947 :     *
948 :     * @return numeric vector of the values in the upper triangles of the
949 :     * Omega matrices
950 :     */
951 :     SEXP ssclme_coef(SEXP x)
952 :     {
953 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
954 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
955 :     i, nf = length(Omega), vind;
956 :     SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
957 :     double *vv = REAL(val);
958 :    
959 :     vind = 0;
960 :     for (i = 0; i < nf; i++) {
961 :     int j, k, nci = nc[i];
962 :     double *omgi = REAL(VECTOR_ELT(Omega, i));
963 :     for (j = 0; j < nci; j++) {
964 :     for (k = 0; k <= j; k++) {
965 :     vv[vind++] = omgi[j*nci + k];
966 :     }
967 :     }
968 :     }
969 :     UNPROTECT(1);
970 :     return val;
971 :     }
972 :    
973 :     /**
974 : bates 19 * Extract the upper triangles of the Omega matrices in the unconstrained
975 :     * parameterization.
976 :     * (These are not in any sense "coefficients" but the extractor is
977 :     * called coef for historical reasons.)
978 :     *
979 :     * @param x pointer to an ssclme object
980 :     *
981 :     * @return numeric vector of the values in the upper triangles of the
982 :     * Omega matrices
983 :     */
984 :     SEXP ssclme_coefUnc(SEXP x)
985 :     {
986 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
987 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
988 :     i, nf = length(Omega), vind;
989 :     SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
990 :     double *vv = REAL(val);
991 :    
992 :     vind = 0;
993 :     for (i = 0; i < nf; i++) {
994 :     int nci = nc[i];
995 :     if (nci == 1) {
996 :     vv[vind++] = log(REAL(VECTOR_ELT(Omega, i))[0]);
997 :     } else {
998 :     int j, k, ncip1 = nci + 1, ncisq = nci * nci;
999 :     double *tmp = Memcpy(Calloc(ncisq, double),
1000 :     REAL(VECTOR_ELT(Omega, i)), ncisq);
1001 :     F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1002 :     if (j) /* should never happen */
1003 :     error("DPOTRF returned error code %d", j);
1004 :     for (j = 0; j < nci; j++) {
1005 :     double diagj = tmp[j * ncip1];
1006 :     vv[vind++] = 2. * log(diagj);
1007 :     for (k = j + 1; k < nci; k++) {
1008 :     tmp[j + k * nci] /= diagj;
1009 :     }
1010 :     }
1011 :     for (j = 0; j < nci; j++) {
1012 :     for (k = j + 1; k < nci; k++) {
1013 :     vv[vind++] = tmp[j + k * nci];
1014 :     }
1015 :     }
1016 :     Free(tmp);
1017 :     }
1018 :     }
1019 :     UNPROTECT(1);
1020 :     return val;
1021 :     }
1022 :    
1023 :     /**
1024 :     * Assign the upper triangles of the Omega matrices in the unconstrained
1025 :     * parameterization.
1026 :     *
1027 :     * @param x pointer to an ssclme object
1028 :     * @param coef pointer to an numeric vector of appropriate length
1029 :     *
1030 :     * @return R_NilValue
1031 :     */
1032 :     SEXP ssclme_coefGetsUnc(SEXP x, SEXP coef)
1033 :     {
1034 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1035 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1036 :     cind, i, nf = length(Omega),
1037 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1038 :     double *cc = REAL(coef);
1039 :    
1040 :     if (length(coef) != coef_length(nf, nc) || !isReal(coef))
1041 :     error("coef must be a numeric vector of length %d",
1042 :     coef_length(nf, nc));
1043 :     cind = 0;
1044 :     for (i = 0; i < nf; i++) {
1045 :     int nci = nc[i];
1046 :     if (nci == 1) {
1047 :     REAL(VECTOR_ELT(Omega, i))[0] = exp(cc[cind++]);
1048 :     } else {
1049 :     int odind = cind + nci, /* off-diagonal index */
1050 :     j, k,
1051 :     ncip1 = nci + 1,
1052 :     ncisq = nci * nci;
1053 :     double
1054 :     *omgi = REAL(VECTOR_ELT(Omega, i)),
1055 :     *tmp = Calloc(ncisq, double),
1056 :     diagj, one = 1.;
1057 :     /* LD in omgi and L' in tmp */
1058 :     memset(omgi, 0, sizeof(double) * ncisq);
1059 :     for (j = 0; j < nci; j++) {
1060 :     omgi[j * ncip1] = diagj = exp(cc[cind++]);
1061 :     for (k = j + 1; k < nci; k++) {
1062 :     omgi[j*nci + k] = diagj * (tmp[k*nci + j] = cc[odind++]);
1063 :     }
1064 :     }
1065 :     F77_CALL(dtrmm)("R", "U", "N", "U", &nci, &nci, &one,
1066 :     tmp, &nci, omgi, &nci);
1067 :     Free(tmp);
1068 :     cind = odind;
1069 :     }
1070 :     }
1071 :     status[0] = status[1] = 0;
1072 :     return x;
1073 :     }
1074 :    
1075 :     /**
1076 : bates 10 * Assign the upper triangles of the Omega matrices.
1077 :     * (These are not in any sense "coefficients" but are
1078 :     * called coef for historical reasons.)
1079 :     *
1080 :     * @param x pointer to an ssclme object
1081 :     * @param coef pointer to an numeric vector of appropriate length
1082 :     *
1083 :     * @return R_NilValue
1084 :     */
1085 :     SEXP ssclme_coefGets(SEXP x, SEXP coef)
1086 :     {
1087 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1088 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1089 :     cind, i, nf = length(Omega),
1090 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1091 :     double *cc = REAL(coef);
1092 :    
1093 :     if (length(coef) != coef_length(nf, nc) || !isReal(coef))
1094 :     error("coef must be a numeric vector of length %d",
1095 :     coef_length(nf, nc));
1096 :     cind = 0;
1097 :     for (i = 0; i < nf; i++) {
1098 :     int j, k, nci = nc[i];
1099 :     double *omgi = REAL(VECTOR_ELT(Omega, i));
1100 :     for (j = 0; j < nci; j++) {
1101 :     for (k = 0; k <= j; k++) {
1102 :     omgi[j*nci + k] = cc[cind++];
1103 :     }
1104 :     }
1105 :     }
1106 :     status[0] = status[1] = 0;
1107 : bates 11 return x;
1108 : bates 10 }
1109 :    
1110 : bates 11 SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1111 : bates 10 {
1112 :     SEXP
1113 :     Omega = GET_SLOT(x, Matrix_OmegaSym),
1114 :     RZXsl = GET_SLOT(x, Matrix_RZXSym),
1115 :     ncsl = GET_SLOT(x, Matrix_ncSym),
1116 :     bVar = GET_SLOT(x, Matrix_bVarSym);
1117 :     int
1118 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1119 :     *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1120 :     *nc = INTEGER(ncsl),
1121 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1122 :     REML = asLogical(REMLp),
1123 :     i, info, iter,
1124 :     n = dims[0],
1125 :     nEM = asInteger(nsteps),
1126 :     nf = length(ncsl) - 2,
1127 :     nobs = nc[nf + 1],
1128 :     p,
1129 :     pp1 = dims[1],
1130 :     verbose = asLogical(verb);
1131 :     double
1132 :     *RZX = REAL(RZXsl),
1133 :     *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),
1134 :     *b,
1135 :     alpha,
1136 :     one = 1.,
1137 :     zero = 0.;
1138 :    
1139 :     p = pp1 - 1;
1140 :     b = RZX + p * n;
1141 :     if (verbose) Rprintf(" EM iterations\n");
1142 :     for (iter = 0; iter <= nEM; iter++) {
1143 :     ssclme_invert(x);
1144 :     if (verbose) {
1145 :     SEXP coef = PROTECT(ssclme_coef(x));
1146 :     int lc = length(coef); double *cc = REAL(coef);
1147 :     Rprintf("%3d %.3f", iter, dev[REML ? 1 : 0]);
1148 :     for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1149 :     Rprintf("\n");
1150 :     UNPROTECT(1);
1151 :     }
1152 :     for (i = 0; i < nf; i++) {
1153 :     int ki = Gp[i+1] - Gp[i],
1154 :     nci = nc[i],
1155 :     mi = ki/nci;
1156 :     double
1157 :     *vali = REAL(VECTOR_ELT(Omega, i));
1158 :    
1159 :     alpha = ((double)(REML?(nobs-p):nobs))/((double)mi);
1160 :     F77_CALL(dsyrk)("U", "N", &nci, &mi,
1161 :     &alpha, b + Gp[i], &nci,
1162 :     &zero, vali, &nci);
1163 :     alpha = 1./((double) mi);
1164 :     F77_CALL(dsyrk)("U", "N", &nci, &ki,
1165 :     &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1166 :     &one, vali, &nci);
1167 :     if (REML) {
1168 :     int mp = mi * p;
1169 :     F77_CALL(dsyrk)("U", "N", &nci, &mp,
1170 :     &alpha, RZX + Gp[i], &nci,
1171 :     &one, vali, &nci);
1172 :     }
1173 :     F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1174 :     if (info) error("DPOTRF returned error code %d", info);
1175 :     F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1176 :     if (info) error("DPOTRF returned error code %d", info);
1177 :     }
1178 :     status[0] = status[1] = 0;
1179 :     }
1180 :     ssclme_factor(x);
1181 :     return R_NilValue;
1182 :     }
1183 : bates 11
1184 :     SEXP ssclme_asSscMatrix(SEXP x)
1185 :     {
1186 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix")));
1187 :     int *dims = INTEGER(GET_SLOT(val, Matrix_DimSym));
1188 :    
1189 :     dims[0] = dims[1] = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
1190 :     SET_SLOT(val, Matrix_pSym, duplicate(GET_SLOT(x, Matrix_pSym)));
1191 :     SET_SLOT(val, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_iSym)));
1192 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
1193 :     CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0))[0] = 'U';
1194 :     UNPROTECT(1);
1195 :     return val;
1196 :     }
1197 :    
1198 :     SEXP ssclme_asTscMatrix(SEXP x)
1199 :     {
1200 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix")));
1201 :     int *dims = INTEGER(GET_SLOT(val, Matrix_DimSym));
1202 :    
1203 :     dims[0] = dims[1] = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
1204 :     SET_SLOT(val, Matrix_pSym, duplicate(GET_SLOT(x, Matrix_LpSym)));
1205 :     SET_SLOT(val, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_LiSym)));
1206 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_LxSym)));
1207 :     CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0))[0] = 'U';
1208 :     CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0] = 'U';
1209 :     UNPROTECT(1);
1210 :     return val;
1211 :     }
1212 :    

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge