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 157 - (view) (download) (as text)

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

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