SCM

SCM Repository

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

Annotation of /pkg/src/lmer.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 415 #include "lmer.h"
2 : bates 429 /* TODO
3 :     * - The egsingle example with ~year|childid+schoolid shows an unusual
4 :     * drop in the deviance when switching from ECME to optim. Is it real?
5 :     * (Apparently so.)
6 :     */
7 : bates 255
8 :     /**
9 : bates 373 * Calculate the length of the parameter vector (historically called "coef"
10 :     * even though these are not coefficients).
11 : bates 255 *
12 :     * @param nf number of factors
13 :     * @param nc number of columns in the model matrices for each factor
14 :     *
15 :     * @return total length of the coefficient vector
16 :     */
17 : bates 298 static R_INLINE
18 : bates 255 int coef_length(int nf, const int nc[])
19 :     {
20 :     int i, ans = 0;
21 :     for (i = 0; i < nf; i++) ans += (nc[i] * (nc[i] + 1))/2;
22 :     return ans;
23 :     }
24 :    
25 : bates 274 /**
26 : bates 466 * Calculate the zero-based index in a row-wise packed lower triangular matrix.
27 :     * This is used for the arrays of blocked sparse matrices.
28 : bates 274 *
29 : bates 275 * @param i column number (zero-based)
30 :     * @param k row number (zero-based)
31 : bates 274 *
32 :     * @return The index of the (k,i) element of a packed lower triangular matrix
33 :     */
34 :     static R_INLINE
35 : bates 447 int Lind(int k, int i)
36 : bates 274 {
37 : bates 447 if (k < i) error("Lind(k = %d, i = %d) must have k >= i", k, i);
38 : bates 466 return (k * (k + 1))/2 + i;
39 : bates 274 }
40 :    
41 :     /**
42 :     * Allocate a 3-dimensional array
43 :     *
44 :     * @param TYP The R Type code (e.g. INTSXP)
45 :     * @param nr number of rows
46 :     * @param nc number of columns
47 :     * @param nf number of faces
48 :     *
49 :     * @return A 3-dimensional array of the indicated dimensions and type
50 :     */
51 : bates 257 static
52 : bates 275 SEXP alloc3Darray(int TYP, int nr, int nc, int nf)
53 : bates 257 {
54 :     SEXP val, dd = PROTECT(allocVector(INTSXP, 3));
55 :    
56 :     INTEGER(dd)[0] = nr; INTEGER(dd)[1] = nc; INTEGER(dd)[2] = nf;
57 :     val = allocArray(TYP, dd);
58 :     UNPROTECT(1);
59 :     return val;
60 :     }
61 :    
62 : bates 298 static R_INLINE
63 :     int match_mat_dims(const int xd[], const int yd[])
64 :     {
65 :     return xd[0] == yd[0] && xd[1] == yd[1];
66 :     }
67 :    
68 : bates 274 /**
69 : bates 415 * Check validity of an lmer object.
70 : bates 274 *
71 : bates 415 * @param x Pointer to an lmer object
72 : bates 274 *
73 : bates 415 * @return TRUE if the object is a valid lmer object, else a string
74 : bates 274 * describing the nature of the violation.
75 :     */
76 : bates 415 SEXP lmer_validate(SEXP x)
77 : bates 255 {
78 : bates 298 SEXP
79 : bates 379 /* ZZxP = GET_SLOT(x, Matrix_ZZxSym), */
80 : bates 298 ZtXP = GET_SLOT(x, Matrix_ZtXSym),
81 :     XtXP = GET_SLOT(x, Matrix_XtXSym),
82 :     RZXP = GET_SLOT(x, Matrix_RZXSym),
83 : bates 387 RXXP = GET_SLOT(x, Matrix_RXXSym)
84 : bates 415 /* , cnames = GET_SLOT(x, Matrix_cnamesSym) */
85 : bates 379 ;
86 : bates 299 int *ZtXd = INTEGER(getAttrib(ZtXP, R_DimSymbol)),
87 :     *XtXd = INTEGER(getAttrib(XtXP, R_DimSymbol));
88 : bates 298
89 :     if (!(isReal(ZtXP) && isReal(XtXP) && isReal(RZXP) && isReal(RXXP) ))
90 :     return ScalarString(mkChar("Slots ZtX, XtX, RZX, and RXX must be real matrices"));
91 :     if (!match_mat_dims(ZtXd, INTEGER(getAttrib(RZXP, R_DimSymbol))))
92 : bates 299 return ScalarString(mkChar("Dimensions of slots ZtX and RZX must match"));
93 :     if (!match_mat_dims(XtXd, INTEGER(getAttrib(RXXP, R_DimSymbol))))
94 :     return ScalarString(mkChar("Dimensions of slots XtX and RXX must match"));
95 :     if (ZtXd[1] != XtXd[0] || XtXd[0] != XtXd[1])
96 :     return ScalarString(mkChar("Slots XtX must be a square matrix with same no. of cols as ZtX"));
97 : bates 255 return ScalarLogical(1);
98 :     }
99 :    
100 : bates 274 static R_INLINE
101 :     int Tind(const int rowind[], const int colptr[], int i, int j)
102 :     {
103 :     int k, k2 = colptr[j + 1];
104 :     for (k = colptr[j]; k < k2; k++)
105 :     if (rowind[k] == i) return k;
106 :     error("row %d and column %d not defined in rowind and colptr",
107 :     i, j);
108 : bates 466 return -1; /* -Wall */
109 : bates 274 }
110 :    
111 : bates 275 /**
112 : bates 415 * Create the pairwise crosstabulation of the elements of flist.
113 : bates 275 *
114 : bates 415 * @param flist pointer to the factor list.
115 :     * @param nobs number of observations.
116 :     * @param nc number of columns in the model matrices.
117 : bates 275 *
118 : bates 415 * @return the pairwise crosstabulation in the form of the ZtZ array.
119 : bates 275 */
120 : bates 415 static SEXP
121 :     lmer_crosstab(SEXP flist, int nobs, const int nc[])
122 : bates 275 {
123 : bates 415 int i, nf = length(flist);
124 :     int npairs = (nf * (nf + 1))/2;
125 :     SEXP val = PROTECT(allocVector(VECSXP, npairs));
126 : bates 485 SEXP cscbCl = PROTECT(MAKE_CLASS("dgBCMatrix"));
127 : bates 415 int *Ti = Calloc(nobs, int),
128 :     *nlevs = Calloc(nf, int),
129 :     **zb = Calloc(nf, int*); /* zero-based indices */
130 :     double *ones = Calloc(nobs, double),
131 :     *Tx = Calloc(nobs, double);
132 : bates 274
133 : bates 415 for (i = 0; i < nobs; i++) ones[i] = 1.;
134 :     for (i = 0; i < nf; i++) { /* populate the zb vectors */
135 :     SEXP fi = VECTOR_ELT(flist, i);
136 :     int j;
137 : bates 274
138 : bates 415 zb[i] = Calloc(nobs, int);
139 :     nlevs[i] = length(getAttrib(fi, R_LevelsSymbol));
140 :     for (j = 0; j < nobs; j++) zb[i][j] = INTEGER(fi)[j] - 1;
141 :     for (j = 0; j <= i; j++) {
142 :     int *ijp, ind = Lind(i, j), nnz;
143 :     SEXP ZZij;
144 :    
145 :     SET_VECTOR_ELT(val, ind, NEW_OBJECT(cscbCl));
146 :     ZZij = VECTOR_ELT(val, ind);
147 :     SET_SLOT(ZZij, Matrix_pSym, allocVector(INTSXP, nlevs[j] + 1));
148 :     ijp = INTEGER(GET_SLOT(ZZij, Matrix_pSym));
149 : bates 488 triplet_to_col(nlevs[i], nlevs[j], nobs, zb[i], zb[j], ones,
150 : bates 415 ijp, Ti, Tx);
151 :     nnz = ijp[nlevs[j]];
152 :     SET_SLOT(ZZij, Matrix_iSym, allocVector(INTSXP, nnz));
153 :     Memcpy(INTEGER(GET_SLOT(ZZij, Matrix_iSym)), Ti, nnz);
154 :     SET_SLOT(ZZij, Matrix_xSym, alloc3Darray(REALSXP, nc[i], nc[j], nnz));
155 :     /* The crosstab counts are copied into the first nnz elements */
156 :     /* of the x slot. These aren't the correct array positions */
157 :     /* unless nc[i] == nc[j] == 1 but we don't use them. */
158 :     Memcpy(REAL(GET_SLOT(ZZij, Matrix_xSym)), Tx, nnz);
159 : bates 274 }
160 :     }
161 : bates 415
162 :     for (i = 0; i < nf; i++) Free(zb[i]);
163 :     Free(zb); Free(nlevs); Free(ones); Free(Ti); Free(Tx);
164 : bates 466 UNPROTECT(2);
165 : bates 255 return val;
166 :     }
167 :    
168 : bates 415
169 : bates 255 /**
170 : bates 415 * Update the arrays ZtZ, ZtX, and XtX in an lme object
171 :     * according to a list of model matrices.
172 : bates 255 *
173 : bates 415 * @param x pointer to an lmer object
174 : bates 255 * @param mmats pointer to a list of model matrices
175 :     *
176 :     * @return NULL
177 :     */
178 : bates 415 SEXP lmer_update_mm(SEXP x, SEXP mmats)
179 : bates 255 {
180 : bates 264 SEXP
181 : bates 415 ZtZP = GET_SLOT(x, Matrix_ZtZSym),
182 : bates 255 ZtXP = GET_SLOT(x, Matrix_ZtXSym),
183 : bates 415 flist = GET_SLOT(x, Matrix_flistSym);
184 :     int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
185 :     *dims = INTEGER(getAttrib(ZtXP, R_DimSymbol)),
186 : bates 296 *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
187 : bates 255 *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
188 : bates 415 nf = length(flist), nfp1 = nf + 1,
189 : bates 255 i, ione = 1,
190 : bates 275 nobs = nc[nfp1],
191 :     pp1 = nc[nf];
192 : bates 255 double
193 :     *X,
194 :     *XtX = REAL(GET_SLOT(x, Matrix_XtXSym)),
195 :     *ZtX = REAL(ZtXP),
196 :     one = 1.0, zero = 0.0;
197 :    
198 : bates 275 if (!isNewList(mmats) || length(mmats) != nfp1)
199 :     error("mmats must be a list of %d model matrices", nfp1);
200 :     for (i = 0; i <= nf; i++) {
201 : bates 255 SEXP mmat = VECTOR_ELT(mmats, i);
202 : bates 415 int *mdims = INTEGER(getAttrib(mmat, R_DimSymbol));
203 : bates 255
204 :     if (!isMatrix(mmat) || !isReal(mmat))
205 :     error("element %d of mmats is not a numeric matrix", i + 1);
206 : bates 415 if (nobs != mdims[0])
207 : bates 255 error("Expected %d rows in the %d'th model matrix. Got %d",
208 : bates 415 nobs, i+1, mdims[0]);
209 :     if (nc[i] != mdims[1])
210 : bates 255 error("Expected %d columns in the %d'th model matrix. Got %d",
211 : bates 415 nc[i], i+1, mdims[1]);
212 : bates 255 }
213 :     /* Create XtX */
214 : bates 275 X = REAL(VECTOR_ELT(mmats, nf));
215 :     F77_CALL(dsyrk)("U", "T", &pp1, &nobs, &one, X, &nobs, &zero, XtX, nc + nf);
216 : bates 255 /* Zero an accumulator */
217 : bates 429 AZERO(ZtX, pp1 * Gp[nf]);
218 : bates 275 for (i = 0; i < nf; i++) {
219 : bates 415 int *fac = INTEGER(VECTOR_ELT(flist, i)),
220 :     j, k, nci = nc[i], ZtXrows = Gp[i+1] - Gp[i];
221 :     int ncisqr = nci * nci, nlev = ZtXrows/nci;
222 :     double *Z = REAL(VECTOR_ELT(mmats, i)), *ZZx;
223 : bates 268
224 :     for (k = 0; k < i; k++) {
225 : bates 415 SEXP ZZxM = VECTOR_ELT(ZtZP, Lind(i, k));
226 : bates 274 int *rowind = INTEGER(GET_SLOT(ZZxM, Matrix_iSym)),
227 :     *colptr = INTEGER(GET_SLOT(ZZxM, Matrix_pSym));
228 : bates 415 int *f2 = INTEGER(VECTOR_ELT(flist, k)), nck = nc[k];
229 : bates 268 double *Zk = REAL(VECTOR_ELT(mmats, k));
230 :    
231 : bates 274 ZZx = REAL(GET_SLOT(ZZxM, Matrix_xSym));
232 : bates 429 AZERO(ZZx, length(GET_SLOT(ZZxM, Matrix_xSym)));
233 : bates 268 for (j = 0; j < nobs; j++) {
234 : bates 274 F77_CALL(dgemm)("T", "N", nc + i, nc + k, &ione, &one,
235 : bates 268 Z + j, &nobs, Zk + j, &nobs, &one,
236 :     ZZx + Tind(rowind, colptr, fac[j] - 1, f2[j] - 1)
237 :     * (nci * nck), &nci);
238 :     }
239 :     }
240 : bates 415 ZZx = REAL(GET_SLOT(VECTOR_ELT(ZtZP, Lind(i, i)), Matrix_xSym));
241 : bates 429 AZERO(ZZx, nci * nci * nlev);
242 : bates 255 if (nci == 1) { /* single column in Z */
243 :     for (j = 0; j < nobs; j++) {
244 :     int fj = fac[j] - 1; /* factor indices are 1-based */
245 :     ZZx[fj] += Z[j] * Z[j];
246 : bates 296 F77_CALL(daxpy)(&pp1, Z + j, X + j, &nobs, ZtX + fj, dims);
247 : bates 255 }
248 :     } else {
249 :     for (j = 0; j < nobs; j++) {
250 :     int fj = fac[j] - 1; /* factor indices are 1-based */
251 :    
252 :     F77_CALL(dsyr)("U", nc + i, &one, Z + j, &nobs,
253 :     ZZx + fj * ncisqr, nc + i);
254 :     F77_CALL(dgemm)("T", "N", nc + i, &pp1, &ione,
255 :     &one, Z + j, &nobs,
256 :     X + j, &nobs, &one,
257 : bates 296 ZtX + fj * nci, dims);
258 : bates 255 }
259 :     }
260 :     ZtX += ZtXrows;
261 :     }
262 :     status[0] = status[1] = 0;
263 :     return R_NilValue;
264 :     }
265 :    
266 : bates 415 /**
267 :     * Create an lmer object from a list grouping factors and a list of model
268 :     * matrices. There is one more model matrix than grouping factor. The last
269 :     * model matrix is the fixed effects and the response.
270 :     *
271 :     * @param facs pointer to a list of grouping factors
272 :     * @param ncv pointer to a list of model matrices
273 :     *
274 :     * @return pointer to an lmer object
275 :     */
276 :     SEXP lmer_create(SEXP flist, SEXP mmats)
277 :     {
278 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("lmer")));
279 :     SEXP ZtZ, cnames, fnms, nms;
280 :     int *nc, i, nf = length(flist), nobs;
281 :    
282 :     /* Check validity of flist */
283 :     if (!(nf > 0 && isNewList(flist)))
284 :     error("flist must be a non-empty list");
285 :     nobs = length(VECTOR_ELT(flist, 0));
286 :     if (nobs < 1) error("flist[[0]] must be a non-null factor");
287 :     for (i = 0; i < nf; i++) {
288 :     SEXP fi = VECTOR_ELT(flist, i);
289 :     if (!(isFactor(fi) && length(fi) == nobs))
290 :     error("flist[[%d]] must be a factor of length %d",
291 :     i + 1, nobs);
292 :     }
293 :     SET_SLOT(val, Matrix_flistSym, flist);
294 :     if (!(isNewList(mmats) && length(mmats) == (nf + 1)))
295 :     error("mmats must be a list of length %d", nf + 1);
296 :     SET_SLOT(val, Matrix_ncSym, allocVector(INTSXP, nf + 2));
297 :     nc = INTEGER(GET_SLOT(val, Matrix_ncSym));
298 :     nc[nf + 1] = nobs;
299 :     for (i = 0; i <= nf; i++) {
300 :     SEXP mi = VECTOR_ELT(mmats, i);
301 :     int *dims;
302 : bates 255
303 : bates 415 if (!(isMatrix(mi) && isReal(mi)))
304 :     error("mmats[[%d]] must be a numeric matrix", i + 1);
305 :     dims = INTEGER(getAttrib(mi, R_DimSymbol));
306 :     if (dims[0] != nobs)
307 :     error("mmats[[%d]] must have %d rows", i + 1, nobs);
308 :     if (dims[1] < 1)
309 :     error("mmats[[%d]] must have at least 1 column", i + 1);
310 :     nc[i] = dims[1];
311 :     } /* Arguments have now been checked for type, dimension, etc. */
312 :     /* Create pairwise crosstabulation in ZtZ */
313 :     SET_SLOT(val, Matrix_ZtZSym, lmer_crosstab(flist, nobs, nc));
314 :     lmer_populate(val);
315 :     ZtZ = GET_SLOT(val, Matrix_ZtZSym);
316 :     /* FIXME: Check for possible reordering of the factors to maximize the
317 : bates 429 * number of levels (columns?) in the leading nested sequence. */
318 : bates 415 fnms = getAttrib(flist, R_NamesSymbol);
319 :     /* Allocate and populate nc and cnames */
320 :     SET_SLOT(val, Matrix_cnamesSym, allocVector(VECSXP, nf + 1));
321 :     cnames = GET_SLOT(val, Matrix_cnamesSym);
322 :     setAttrib(cnames, R_NamesSymbol, allocVector(STRSXP, nf + 1));
323 :     nms = getAttrib(cnames, R_NamesSymbol);
324 :     for (i = 0; i <= nf; i++) {
325 :     SEXP mi = VECTOR_ELT(mmats, i);
326 :     SET_VECTOR_ELT(cnames, i,
327 :     duplicate(VECTOR_ELT(getAttrib(mi, R_DimNamesSymbol),
328 :     1)));
329 :     if (i < nf)
330 :     SET_STRING_ELT(nms, i, duplicate(STRING_ELT(fnms, i)));
331 :     else
332 :     SET_STRING_ELT(nms, nf, mkChar(".fixed"));
333 :     }
334 :     lmer_update_mm(val, mmats);
335 : bates 429 SET_SLOT(val, Matrix_bVarSym, duplicate(GET_SLOT(val, Matrix_DSym)));
336 : bates 415 UNPROTECT(1);
337 :     return val;
338 :     }
339 :    
340 : bates 255 /**
341 :     * Create and insert initial values for Omega.
342 :     *
343 : bates 415 * @param x pointer to an lmer object
344 : bates 255 *
345 :     * @return NULL
346 :     */
347 : bates 415 SEXP lmer_initial(SEXP x)
348 : bates 255 {
349 : bates 415 SEXP Omg = GET_SLOT(x, Matrix_OmegaSym);
350 :     int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)), i, nf = length(Omg);
351 : bates 268
352 :     for (i = 0; i < nf; i++) {
353 : bates 415 SEXP ZZxP = GET_SLOT(VECTOR_ELT(GET_SLOT(x, Matrix_ZtZSym), Lind(i, i)),
354 : bates 275 Matrix_xSym);
355 : bates 415 int *dims = INTEGER(getAttrib(ZZxP, R_DimSymbol));
356 :     int j, k, nzc = dims[0], nlev = dims[2];
357 :     int nzcsqr = nzc * nzc, nzcp1 = nzc + 1;
358 :     double *Omega = REAL(VECTOR_ELT(Omg, i)),
359 : bates 275 mi = 0.375 / ((double) nlev);
360 : bates 255
361 : bates 429 AZERO(Omega, nzc * nzc);
362 : bates 268 for (j = 0; j < nlev; j ++) {
363 :     for (k = 0; k < nzc; k++) {
364 : bates 275 Omega[k * nzcp1] += REAL(ZZxP)[k * nzcp1 + j * nzcsqr] * mi;
365 : bates 268 }
366 : bates 255 }
367 :     }
368 :     status[0] = status[1] = 0;
369 :     return R_NilValue;
370 :     }
371 :    
372 :     /**
373 : bates 437 * Copy ZtZ to ZZpO and L. Inflate diagonal blocks of ZZpO by Omega.
374 : bates 415 * Update devComp[1].
375 : bates 300 *
376 : bates 415 * @param x pointer to an lmer object
377 : bates 300 */
378 : bates 415 SEXP
379 :     lmer_inflate(SEXP x)
380 : bates 300 {
381 : bates 415 SEXP Omg = GET_SLOT(x, Matrix_OmegaSym),
382 :     ZZpO = GET_SLOT(x, Matrix_ZZpOSym),
383 : bates 445 ZtZ = GET_SLOT(x, Matrix_ZtZSym),
384 :     LP = GET_SLOT(x, Matrix_LSym);
385 : bates 415 int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
386 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
387 :     i, k, nf = length(Omg);
388 :     double *dcmp = REAL(GET_SLOT(x, Matrix_devCompSym));
389 : bates 300
390 : bates 415 for (i = 0; i < nf; i++) {
391 : bates 445 SEXP ZZOel = VECTOR_ELT(ZZpO, i);
392 :     SEXP ZZOm = GET_SLOT(ZZOel, Matrix_xSym);
393 :     SEXP ZZel = VECTOR_ELT(ZtZ, Lind(i, i));
394 :     int *Di = INTEGER(GET_SLOT(ZZOel, Matrix_iSym)),
395 :     *Dp = INTEGER(GET_SLOT(ZZOel, Matrix_pSym)),
396 :     *Si = INTEGER(GET_SLOT(ZZel, Matrix_iSym)),
397 :     *Sp = INTEGER(GET_SLOT(ZZel, Matrix_pSym)),
398 :     *dims = INTEGER(getAttrib(ZZOm, R_DimSymbol));
399 :     int sz = dims[0] * dims[1];
400 :     int ii, j, nci = nc[i], ncisqr = nci * nci;
401 : bates 415 int nlev = (Gp[i + 1] - Gp[i])/nci;
402 : bates 445 double *Omega = REAL(VECTOR_ELT(Omg, i)),
403 :     *ZZ = REAL(GET_SLOT(ZZel, Matrix_xSym)),
404 :     *tmp = Memcpy(Calloc(ncisqr, double), Omega, ncisqr);
405 : bates 300
406 : bates 415 F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j); /* update dcmp[1] */
407 :     if (j)
408 :     error("Leading %d minor of Omega[[%d]] not positive definite",
409 :     j, i + 1);
410 :     for (j = 0; j < nci; j++) { /* nlev * logDet(Omega_i) */
411 :     dcmp[1] += nlev * 2. * log(tmp[j * (nci + 1)]);
412 : bates 300 }
413 : bates 415 Free(tmp);
414 : bates 447 AZERO(REAL(ZZOm), dims[0] * dims[1] * dims[2]);
415 : bates 445 for (j = 0; j < nlev; j++) { /* copy diagonal block and inflate */
416 :     double *ZZOkk = REAL(ZZOm) + Tind(Di, Dp, j, j) * sz;
417 :     int kk, k2 = Sp[j + 1];
418 :     for (kk = Sp[j]; kk < k2; kk++) {
419 :     Memcpy(REAL(ZZOm) + Tind(Di, Dp, Si[kk], j) * sz,
420 :     ZZ + kk * sz, sz);
421 :     }
422 :     for (kk = 0; kk < nci; kk++) {
423 :     for (ii = 0; ii <= kk; ii++) {
424 :     int ind = ii + kk * nci;
425 :     ZZOkk[ind] += Omega[ind];
426 :     }
427 :     }
428 :     }
429 :     for (k = i + 1; k < nf; k++) {
430 : bates 415 int ind = Lind(k, i);
431 : bates 445 SEXP Lel = VECTOR_ELT(LP, ind),
432 :     Lm = GET_SLOT(Lel, Matrix_xSym);
433 :     double *L = REAL(Lm);
434 : bates 300
435 : bates 445 dims = INTEGER(getAttrib(Lm, R_DimSymbol));
436 :     ZZel = VECTOR_ELT(ZtZ, ind);
437 :     ZZ = REAL(GET_SLOT(ZZel, Matrix_xSym));
438 :     Di = INTEGER(GET_SLOT(Lel, Matrix_iSym));
439 :     Dp = INTEGER(GET_SLOT(Lel, Matrix_pSym));
440 :     Si = INTEGER(GET_SLOT(ZZel, Matrix_iSym));
441 :     Sp = INTEGER(GET_SLOT(ZZel, Matrix_pSym));
442 :     sz = dims[0] * dims[1];
443 :    
444 :     AZERO(L, sz * dims[2]); /* zero L */
445 : bates 429 for (j = 0; j < nlev; j++) { /* copy src blocks to dest */
446 : bates 415 int kk, k2 = Sp[j + 1];
447 :     for (kk = Sp[j]; kk < k2; kk++) {
448 : bates 445 Memcpy(L + Tind(Di, Dp, Si[kk], j) * sz, ZZ + kk * sz, sz);
449 : bates 415 }
450 : bates 300 }
451 : bates 274 }
452 :     }
453 : bates 415 return R_NilValue;
454 : bates 274 }
455 :    
456 : bates 445 /**
457 :     * Convert the extended parent pair (Parent, Block) to a parent array
458 :     * for the jth diagonal block of size n.
459 :     *
460 :     * @param j index (0-based) of the diagonal outer block
461 :     * @param n number of inner column blocks in the outer block
462 :     * @param par array of length n to be filled with the parent array
463 :     * @param Parent parent index in the extended parent pair
464 :     * @param Block block index in the extended parent pair
465 :     *
466 :     * @return par
467 :     */
468 : bates 447 static R_INLINE
469 : bates 445 int *block_parent(int j, int n, int par[], const int Parent[], const int Block[])
470 :     {
471 :     int i;
472 :     for (i = 0; i < n; i++) par[i] = (Block[i] == j) ? Parent[i] : -1;
473 :     return par;
474 :     }
475 : bates 291
476 :     /**
477 : bates 415 * If status[["factored"]] is FALSE, create and factor Z'Z+Omega. Also
478 : bates 255 * create RZX and RXX, the deviance components, and the value of the
479 :     * deviance for both ML and REML.
480 :     *
481 : bates 415 * @param x pointer to an lmer object
482 : bates 255 *
483 :     * @return NULL
484 :     */
485 : bates 415 SEXP lmer_factor(SEXP x)
486 : bates 255 {
487 :     int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
488 :    
489 :     if (!status[0]) {
490 : bates 415 SEXP DP = GET_SLOT(x, Matrix_DSym),
491 : bates 274 LP = GET_SLOT(x, Matrix_LSym),
492 : bates 437 RZXP = GET_SLOT(x, Matrix_RZXSym),
493 : bates 415 ZZOP = GET_SLOT(x, Matrix_ZZpOSym),
494 :     Parent = GET_SLOT(x, Matrix_ParentSym);
495 : bates 437 int *dims = INTEGER(getAttrib(RZXP, R_DimSymbol)),
496 : bates 255 *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
497 : bates 415 *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
498 :     i, j, nf = length(DP);
499 : bates 255 int nml = nc[nf + 1], nreml = nml + 1 - nc[nf];
500 :     double
501 :     *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
502 : bates 437 *RZX = REAL(RZXP),
503 : bates 255 *dcmp = REAL(GET_SLOT(x, Matrix_devCompSym)),
504 :     *deviance = REAL(GET_SLOT(x, Matrix_devianceSym)),
505 :     minus1 = -1., one = 1.;
506 :    
507 :    
508 : bates 415 dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;
509 :     Memcpy(RZX, REAL(GET_SLOT(x, Matrix_ZtXSym)), dims[0] * dims[1]);
510 :     Memcpy(RXX, REAL(GET_SLOT(x, Matrix_XtXSym)), dims[1] * dims[1]);
511 : bates 445 lmer_inflate(x); /* initialize ZZpO and L */
512 : bates 255 for (i = 0; i < nf; i++) {
513 : bates 445 SEXP ZZOiP = VECTOR_ELT(ZZOP, i);
514 : bates 415 SEXP DiP = VECTOR_ELT(DP, i);
515 : bates 447 SEXP LiP = VECTOR_ELT(LP, Lind(i, i));
516 : bates 415 int nlev = INTEGER(getAttrib(DiP, R_DimSymbol))[2];
517 :     int jj, nci = nc[i], ncisqr = nci * nci;
518 : bates 445 int *Pari = block_parent(i, nlev, Calloc(nlev, int),
519 :     INTEGER(VECTOR_ELT(VECTOR_ELT(Parent, i), 0)),
520 :     INTEGER(VECTOR_ELT(VECTOR_ELT(Parent, i), 1)));
521 : bates 415 double *D = REAL(DiP);
522 : bates 255
523 : bates 437 jj = cscb_ldl(ZZOiP, Pari, LiP, DiP);
524 :     if (jj != nlev) error("cscb_ldl returned %d < nlev = %d", jj, nlev);
525 : bates 429 for (j = 0; j < nlev; j++) { /* accumulate dcmp[0] */
526 : bates 415 double *Dj = D + j * ncisqr;
527 :     for (jj = 0; jj < nci; jj++) /* accumulate determinant */
528 :     dcmp[0] += 2. * log(Dj[jj * (nci + 1)]);
529 :     }
530 : bates 429 /* Solve L_{i,i} %*% RZX_i := RZX_i */
531 : bates 447 cscb_trsm(LFT, NTR, UNT, 1., LiP, RZX + Gp[i],
532 : bates 415 Gp[i+1] - Gp[i], dims[1], dims[0]);
533 : bates 445 /* Solve D_i^{T/2} %*% RZX_i := RZX_i */
534 : bates 415 for (jj = 0; jj < nlev; jj++) {
535 :     F77_CALL(dtrsm)("L", "U", "T", "N", &nci, dims + 1,
536 :     &one, D + jj * ncisqr, &nci,
537 :     RZX + Gp[i] + jj * nci, dims);
538 :     }
539 :     for (j = i + 1; j < nf; j++) { /* further blocks */
540 : bates 445 SEXP Lji = VECTOR_ELT(LP, Lind(j, i));
541 :     SEXP Lx = GET_SLOT(Lji, Matrix_xSym);
542 :     double *L = REAL(Lx);
543 :     int *xdims = INTEGER(getAttrib(Lx, R_DimSymbol)),
544 :     *Lp = INTEGER(GET_SLOT(Lji, Matrix_pSym));
545 : bates 415 int ntot = xdims[0] * xdims[1];
546 :    
547 : bates 445 /* L_{j,i} := L_{j,i} %*% L_{i,i}^{-T} %*% D_i^{-1/2} */
548 : bates 447 cscb_trcbsm(RGT, LOW, TRN, UNT, 1.0, LiP, Pari, Lji);
549 : bates 415 for (jj = 0; jj < nlev; jj++) {
550 : bates 445 int k, k2 = Lp[jj + 1];
551 :     for (k = Lp[jj]; k < k2; k++)
552 : bates 422 F77_CALL(dtrsm)("R", "U", "N", "N", xdims, xdims + 1,
553 : bates 415 &one, D + jj * ncisqr, &nci,
554 : bates 445 L + k * ntot, xdims);
555 : bates 373 }
556 : bates 445 /* RZX_j := RZX_j - L_{j,i} %*% RZX_i */
557 : bates 447 cscb_mm(LFT, NTR, Gp[j + 1] - Gp[j], dims[1], Gp[i+1] - Gp[i],
558 : bates 445 -1.0, Lji, RZX + Gp[i], dims[0],
559 : bates 415 1.0, RZX + Gp[j], dims[0]);
560 :     }
561 :     for (j = i + 1; j < nf; j++) { /* block pairs and final update */
562 : bates 445 SEXP Lji = VECTOR_ELT(LP, Lind(j, i));
563 :     SEXP Lx = GET_SLOT(Lji, Matrix_xSym);
564 :     double *L = REAL(Lx);
565 :     int *xdims = INTEGER(getAttrib(Lx, R_DimSymbol)),
566 :     *Lp = INTEGER(GET_SLOT(Lji, Matrix_pSym));
567 : bates 415 int ntot = xdims[0] * xdims[1];
568 :    
569 :    
570 : bates 445 /* ZZpO_{j,j} := ZZpO_{j,j} - L{j,i} %*% L_{j,i}^T */
571 : bates 447 cscb_syrk(UPP, NTR, -1.0, Lji, 1.0, VECTOR_ELT(ZZOP, j));
572 : bates 415 for (jj = j+1; jj < nf; jj++) {
573 : bates 445 /* L_{jj,j} := L_{jj,j} - L{jj,i} %*% L_{j,i}^T */
574 : bates 447 cscb_cscbm(NTR, TRN, -1.0, VECTOR_ELT(LP, Lind(jj, i)),
575 : bates 445 Lji, 1.0, VECTOR_ELT(LP, Lind(jj, j)));
576 : bates 255 }
577 : bates 445 /* L_{j,i} := L_{j,i} %*% D_i^{-T/2} */
578 : bates 415 for (jj = 0; jj < nlev; jj++) {
579 : bates 445 int k, k2 = Lp[jj + 1];
580 :     for (k = Lp[jj]; k < k2; k++)
581 : bates 422 F77_CALL(dtrsm)("R", "U", "T", "N", xdims, xdims + 1,
582 : bates 415 &one, D + jj * ncisqr, &nci,
583 : bates 445 L + k * ntot, xdims);
584 : bates 415 }
585 : bates 255 }
586 : bates 445 Free(Pari);
587 : bates 301 }
588 :     /* downdate and factor XtX */
589 : bates 255 F77_CALL(dsyrk)("U", "T", dims + 1, dims,
590 : bates 415 &minus1, RZX, dims, &one, RXX, dims + 1);
591 : bates 255 F77_CALL(dpotrf)("U", dims + 1, RXX, dims + 1, &j);
592 :     if (j) {
593 :     warning("Leading minor of size %d of downdated X'X is indefinite",
594 : bates 373 j);
595 : bates 255 dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;
596 :     } else {
597 :     for (j = 0; j < (dims[1] - 1); j++) /* 2 logDet(RXX) */
598 :     dcmp[2] += 2 * log(RXX[j * (dims[1] + 1)]);
599 :     dcmp[3] = 2. * log(RXX[dims[1] * dims[1] - 1]); /* 2 log(ryy) */
600 :     deviance[0] = /* ML criterion */
601 :     dcmp[0] - dcmp[1] + nml*(1.+dcmp[3]+log(2.*PI/nml));
602 :     deviance[1] = dcmp[0] - dcmp[1] + /* REML */
603 :     dcmp[2] + nreml*(1.+dcmp[3]+log(2.*PI/nreml));
604 :     }
605 :     status[0] = 1; status[1] = 0; /* factored but not inverted */
606 :     }
607 :     return R_NilValue;
608 :     }
609 :    
610 :     /**
611 : bates 415 * Solve one of the matrix equations op(L)*X=alpha*B or
612 :     * X*op(L)=alpha*B where L is a sparse, blocked, unit lower triangular matrix.
613 :     *
614 :     * @param side 'L' for left, 'R' for right
615 :     * @param trans 'T' for transpose, otherwise no transpose
616 :     * @param nf number of grouping factors
617 :     * @param Gp group pointers for the rows
618 :     * @param n number of columns
619 :     * @param alpha multiplier
620 : bates 437 * @param L pointer to the L cscb object
621 : bates 415 * @param mm pointer to the matrix of right-hand sides
622 :     */
623 :     static void
624 :     lmer_sm(char side, char trans, int nf, const int Gp[], int n,
625 : bates 445 double alpha, SEXP L, double B[], int ldb)
626 : bates 415 {
627 :     int itr = (trans == 'T' || trans == 't'), j, k,
628 :     lside = (side == 'L' || side == 'l');
629 :    
630 :     if (lside) {
631 :     if (itr) {
632 :     for (j = nf - 1; j >= 0; j--) {
633 :     int nrj = Gp[j + 1] - Gp[j];
634 :    
635 : bates 447 cscb_trsm(LOW, TRN, UNT, alpha, VECTOR_ELT(L, Lind(j, j)),
636 : bates 415 B + Gp[j], nrj, n, ldb);
637 :     for (k = 0; k < j; k++) {
638 : bates 447 cscb_mm(LFT, TRN, Gp[k + 1] - Gp[k], n, nrj,
639 : bates 445 -1., VECTOR_ELT(L, Lind(j, k)),
640 : bates 415 B + Gp[j], ldb, alpha, B + Gp[k], ldb);
641 :     }
642 :     }
643 :     } else error("Code for non-transpose case not yet written");
644 :     } else error("Code for right-side solutions not yet written");
645 :     }
646 :    
647 :     /**
648 : bates 443 * Fill the nnz array with the number of nonzero inner blocks in each
649 : bates 447 * outer block of the jth inner column block of the ith outer block of
650 :     * L^{-1}. Also allocate the tmp and ind arrays and fill the ind array.
651 : bates 443 *
652 :     * @param i outer block index
653 :     * @param j inner block index within the ith outer block
654 :     * @param nf number of factors
655 :     * @param Parent pointer to the extended parent pairs
656 : bates 447 * @param nc
657 : bates 443 * @param nnz array of length nf
658 : bates 447 * @param tmp array of length nf of pointers to doubles
659 :     * @param ind array of length nf of pointers to ints
660 : bates 443 *
661 :     */
662 : bates 447 #define BLK(i,j) INTEGER(VECTOR_ELT(VECTOR_ELT(Parent, i), 1))[j]
663 :     #define PAR(i,j) INTEGER(VECTOR_ELT(VECTOR_ELT(Parent, i), 0))[j]
664 :    
665 : bates 443 static
666 : bates 447 void fill_nnz(int i, int j, int nf, SEXP Parent, const int nc[],
667 :     int nnz[], double *tmp[], int *ind[])
668 : bates 443 {
669 : bates 447 int *ik = Calloc(nf, int), blk, k, par;
670 : bates 443
671 : bates 447 AZERO(nnz, nf);
672 :     for (blk = BLK(i,j), par = PAR(i,j); blk >= 0;
673 :     k = BLK(blk,par), par = PAR(blk,par), blk = k) nnz[blk]++;
674 :     for (k = 0; k < nf; k++) {
675 :     if (nnz[k]) {
676 :     int sz = nc[i] * nc[k];
677 :     tmp[k] = Calloc(sz * nnz[k], double);
678 :     AZERO(tmp[k], sz * nnz[k]);
679 :     ind[k] = Calloc(nnz[k], int);
680 :     } else {
681 :     tmp[k] = (double *) NULL;
682 :     ind[k] = (int *) NULL;
683 :     }
684 : bates 443 }
685 : bates 447 AZERO(ik, nf);
686 :     for (blk = BLK(i,j), par = PAR(i,j); blk >= 0;
687 :     k = BLK(blk,par), par = PAR(blk,par), blk = k) {
688 :     ind[blk][ik[blk]++] = par;
689 :     }
690 :     Free(ik);
691 : bates 443 }
692 : bates 447
693 :     static R_INLINE
694 :     int fsrch(int target, const int vals[], int nvals)
695 :     {
696 :     int i;
697 :     for (i = 0; i < nvals; i++) if (vals[i] == target) return i;
698 :     error("fsrch: unable to find target %d in nvals %d ", target, nvals);
699 :     return -1; /* -Wall */
700 :     }
701 :    
702 : bates 443 /**
703 : bates 255 * If necessary, factor Z'Z+Omega, ZtX, and XtX then, if necessary,
704 :     * replace the RZX and RXX slots by the corresponding parts of the
705 :     * inverse of the Cholesky factor. Replace the elements of the D slot
706 : bates 429 * by the blockwise inverses and evaluate bVar.
707 : bates 255 *
708 : bates 415 * @param x pointer to an lmer object
709 : bates 255 *
710 :     * @return NULL (x is updated in place)
711 :     */
712 : bates 415 SEXP lmer_invert(SEXP x)
713 : bates 255 {
714 :     int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
715 : bates 415 if (!status[0]) lmer_factor(x);
716 : bates 255 if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))
717 :     error("Unable to invert singular factor of downdated X'X");
718 :     if (!status[1]) {
719 : bates 429 SEXP DP = GET_SLOT(x, Matrix_DSym),
720 :     LP = GET_SLOT(x, Matrix_LSym),
721 : bates 437 ParP = GET_SLOT(x, Matrix_ParentSym),
722 : bates 429 RZXP = GET_SLOT(x, Matrix_RZXSym),
723 :     bVarP = GET_SLOT(x, Matrix_bVarSym);
724 : bates 415 int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
725 : bates 429 *dims = INTEGER(getAttrib(RZXP, R_DimSymbol)),
726 : bates 255 *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
727 : bates 445 i, nf = length(DP);
728 : bates 447 int **ind = Calloc(nf, int *),
729 : bates 443 *nnz = Calloc(nf, int);
730 : bates 447 double **tmp = Calloc(nf, double *),
731 : bates 443 *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
732 : bates 447 *RZX = REAL(RZXP),
733 : bates 443 minus1 = -1., one = 1., zero = 0.;
734 :    
735 : bates 422 /* RXX := RXX^{-1} */
736 : bates 255 F77_CALL(dtrtri)("U", "N", dims + 1, RXX, dims + 1, &i);
737 :     if (i)
738 :     error("Leading minor of size %d of downdated X'X,is indefinite",
739 :     i + 1);
740 : bates 443
741 : bates 422 /* RZX := - RZX %*% RXX */
742 : bates 255 F77_CALL(dtrmm)("R", "U", "N", "N", dims, dims + 1, &minus1,
743 :     RXX, dims + 1, RZX, dims);
744 :     for(i = 0; i < nf; i++) {
745 : bates 447 int info, j, jj, nci = nc[i];
746 : bates 415 int ncisqr = nci * nci, nlev = (Gp[i+1] - Gp[i])/nci;
747 : bates 429 double *Di = REAL(VECTOR_ELT(DP, i)),
748 : bates 415 *RZXi = RZX + Gp[i];
749 : bates 255
750 : bates 422 /* D_i := D_i^{-1}; RZX_i := D_i %*% RZX_i */
751 : bates 255 if (nci == 1) {
752 :     for (j = 0; j < nlev; j++) {
753 :     Di[j] = 1./Di[j];
754 :     for (jj = 0; jj < dims[1]; jj++)
755 : bates 415 RZXi[j + jj * dims[0]] *= Di[j];
756 : bates 255 }
757 :     } else {
758 :     for (j = 0; j < nlev; j++) {
759 :     F77_CALL(dtrtri)("U", "N", &nci, Di + j * ncisqr, &nci, &info);
760 :     if (info)
761 :     error("D[,,%d] for factor %d is singular", j + 1, i + 1);
762 :     F77_CALL(dtrmm)("L", "U", "N", "N", &nci, dims + 1, &one,
763 : bates 415 Di + j * ncisqr, &nci, RZXi + j * nci, dims);
764 : bates 255 }
765 :     }
766 :     }
767 : bates 443
768 : bates 422 /* RZX := L^{-T} %*% RZX */
769 : bates 445 lmer_sm('L', 'T', nf, Gp, dims[1], 1.0, LP, RZX, dims[0]);
770 : bates 437
771 : bates 443 /* Create bVar arrays as crossprod of column blocks of D^{-T/2}%*%L^{-1} */
772 : bates 429 for (i = 0; i < nf; i++) {
773 : bates 447 int j, k, kj, nci = nc[i];
774 : bates 429 int ncisqr = nci * nci, nlev = (Gp[i+1] - Gp[i])/nci;
775 : bates 443 double *Di = REAL(VECTOR_ELT(DP, i)), *bVi = REAL(VECTOR_ELT(bVarP, i));
776 :    
777 : bates 447 AZERO(bVi, nlev * ncisqr); /* zero the accumulator */
778 : bates 429 for (j = 0; j < nlev; j++) {
779 : bates 437 double *bVij = bVi + j * ncisqr, *Dij = Di + j * ncisqr;
780 : bates 443
781 : bates 437 F77_CALL(dsyrk)("U", "N", &nci, &nci, &one, Dij, &nci, &zero, bVij, &nci);
782 : bates 447 /* count non-zero blocks; allocate and zero storage */
783 :     fill_nnz(i, j, nf, ParP, nc, nnz, tmp, ind);
784 :     /* kth row of outer blocks */
785 :     for (k = i; k < nf; k++) {
786 :     SEXP Lki = VECTOR_ELT(LP, Lind(k, i));
787 :     int *Lkii = INTEGER(GET_SLOT(Lki, Matrix_iSym)),
788 :     *Lkip = INTEGER(GET_SLOT(Lki, Matrix_pSym));
789 :     double *Lkix = REAL(GET_SLOT(Lki, Matrix_xSym));
790 :     int kk, sz = nc[i] * nc[k];
791 :    
792 :     /* initialize tmp from jth column of (k,i)th block */
793 :     /* - sign in sol'n incorporated in dtrmm call below */
794 :     for (kk = Lkip[j]; kk < Lkip[j + 1]; kk++)
795 :     Memcpy(tmp[k] + fsrch(Lkii[kk], ind[k], nnz[k]) * sz,
796 :     Lkix + kk * sz, sz);
797 :     /* columns in ind[kk] for (k,kk)th block */
798 :     for (kk = i; kk <= k; kk++) {
799 :     int szk = nc[k] * nc[kk];
800 : bates 445
801 : bates 447 if (!nnz[kk]) continue; /* skip getting slots if not using them */
802 :     Lki = VECTOR_ELT(LP, Lind(k, kk));
803 :     Lkii = INTEGER(GET_SLOT(Lki, Matrix_iSym));
804 :     Lkip = INTEGER(GET_SLOT(Lki, Matrix_pSym));
805 :     Lkix = REAL(GET_SLOT(Lki, Matrix_xSym));
806 :     for (kj = 0; kj < nnz[kk]; kj++) {
807 :     int col = ind[kk][kj], k1;
808 :    
809 :     for (k1 = Lkip[col]; k1 < Lkip[col + 1]; k1++) {
810 :     if ((kk == k) && col >= Lkii[k1]) break;
811 :     F77_CALL(dgemm)("N", "N", nc+k, &nci, nc+kk,
812 :     &minus1, Lkix + k1 * szk, nc + k,
813 :     tmp[kk] + kj * szk, nc+k, &one,
814 :     tmp[k]+fsrch(Lkii[k1],ind[k],nnz[k])*sz,
815 :     nc + k);
816 :     }
817 : bates 437 }
818 : bates 429 }
819 : bates 443 }
820 :     for (k = 0; k < nf; k++) {
821 : bates 447 for (kj = 0; kj < nnz[k]; kj++) {
822 :     F77_CALL(dtrmm)("L", "U", "T", "N", nc + k, &nci, &minus1,
823 :     REAL(VECTOR_ELT(DP, k))+ind[k][kj]*nc[k]*nc[k],
824 :     nc + k, tmp[k] + kj * nc[i] * nc[k], nc + k);
825 : bates 437 }
826 : bates 443 if (nnz[k] > 0) {
827 : bates 447 kj = nc[k] * nnz[k];
828 :     F77_CALL(dsyrk)("U", "T", &nci, &kj, &one, tmp[k], &kj,
829 : bates 443 &one, bVij, &nci);
830 :     }
831 : bates 429 }
832 :     }
833 : bates 443 for (k = 0; k < nf; k++) {
834 :     if (tmp[k]) Free(tmp[k]);
835 :     if (ind[k]) Free(ind[k]);
836 :     }
837 : bates 429 }
838 : bates 447 Free(tmp); Free(nnz); Free(ind);
839 : bates 255 }
840 : bates 443 status[1] = 1;
841 : bates 255 return R_NilValue;
842 :     }
843 :    
844 :     /**
845 :     * Extract the ML or REML conditional estimate of sigma
846 :     *
847 :     * @param x pointer to an lme object
848 :     * @param REML logical scalar - TRUE if REML estimates are requested
849 :     *
850 :     * @return pointer to a numeric scalar
851 :     */
852 : bates 415 SEXP lmer_sigma(SEXP x, SEXP REML)
853 : bates 255 {
854 :     SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
855 :     int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],
856 :     nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))
857 :     [length(GET_SLOT(x, Matrix_OmegaSym)) + 1];
858 :    
859 : bates 415 lmer_invert(x);
860 : bates 255 return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *
861 :     sqrt((double)(asLogical(REML) ?
862 :     nobs + 1 - pp1 : nobs))));
863 :     }
864 :    
865 :     /**
866 :     * Extract the upper triangles of the Omega matrices. These aren't
867 :     * "coefficients" but the extractor is called coef for historical
868 :     * reasons. Within each group these values are in the order of the
869 :     * diagonal entries first then the strict upper triangle in row
870 :     * order.
871 :     *
872 :     * @param x pointer to an lme object
873 :     * @param Unc pointer to a logical scalar indicating if the parameters
874 :     * are in the unconstrained form.
875 :     *
876 :     * @return numeric vector of the values in the upper triangles of the
877 :     * Omega matrices
878 :     */
879 : bates 415 SEXP lmer_coef(SEXP x, SEXP Unc)
880 : bates 255 {
881 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
882 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
883 :     i, nf = length(Omega), unc = asLogical(Unc), vind;
884 :     SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
885 :     double *vv = REAL(val);
886 :    
887 :     vind = 0; /* index in vv */
888 :     for (i = 0; i < nf; i++) {
889 :     int nci = nc[i], ncip1 = nci + 1;
890 :     if (nci == 1) {
891 :     vv[vind++] = (unc ?
892 :     log(REAL(VECTOR_ELT(Omega, i))[0]) :
893 :     REAL(VECTOR_ELT(Omega, i))[0]);
894 :     } else {
895 :     if (unc) { /* L log(D) L' factor of Omega[,,i] */
896 :     int j, k, ncisq = nci * nci;
897 :     double *tmp = Memcpy(Calloc(ncisq, double),
898 :     REAL(VECTOR_ELT(Omega, i)), ncisq);
899 :     F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
900 :     if (j) /* should never happen */
901 :     error("DPOTRF returned error code %d on Omega[[%d]]",
902 :     j, i+1);
903 :     for (j = 0; j < nci; j++) {
904 :     double diagj = tmp[j * ncip1];
905 :     vv[vind++] = 2. * log(diagj);
906 :     for (k = j + 1; k < nci; k++) {
907 :     tmp[j + k * nci] /= diagj;
908 :     }
909 :     }
910 :     for (j = 0; j < nci; j++) {
911 :     for (k = j + 1; k < nci; k++) {
912 :     vv[vind++] = tmp[j + k * nci];
913 :     }
914 :     }
915 :     Free(tmp);
916 :     } else { /* upper triangle of Omega[,,i] */
917 :     int j, k, odind = vind + nci;
918 :     double *omgi = REAL(VECTOR_ELT(Omega, i));
919 :    
920 :     for (j = 0; j < nci; j++) {
921 :     vv[vind++] = omgi[j * ncip1];
922 :     for (k = j + 1; k < nci; k++) {
923 :     vv[odind++] = omgi[k*nci + j];
924 :     }
925 :     }
926 :     vind = odind;
927 :     }
928 :     }
929 :     }
930 :     UNPROTECT(1);
931 :     return val;
932 :     }
933 :    
934 :    
935 :     /**
936 :     * Assign the upper triangles of the Omega matrices.
937 :     * (Called coef for historical reasons.)
938 :     *
939 :     * @param x pointer to an lme object
940 :     * @param coef pointer to an numeric vector of appropriate length
941 :     * @param Unc pointer to a logical scalar indicating if the parameters
942 :     * are in the unconstrained form.
943 :     *
944 :     * @return R_NilValue
945 :     */
946 : bates 415 SEXP lmer_coefGets(SEXP x, SEXP coef, SEXP Unc)
947 : bates 255 {
948 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
949 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
950 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
951 :     cind, i, nf = length(Omega),
952 :     unc = asLogical(Unc);
953 :     double *cc = REAL(coef);
954 :    
955 :     if (length(coef) != coef_length(nf, nc) || !isReal(coef))
956 :     error("coef must be a numeric vector of length %d",
957 :     coef_length(nf, nc));
958 :     cind = 0;
959 :     for (i = 0; i < nf; i++) {
960 :     int nci = nc[i];
961 :     if (nci == 1) {
962 :     REAL(VECTOR_ELT(Omega, i))[0] = (unc ?
963 :     exp(cc[cind++]) :
964 :     cc[cind++]);
965 :     } else {
966 :     int odind = cind + nci, /* off-diagonal index */
967 :     j, k,
968 :     ncip1 = nci + 1,
969 :     ncisq = nci * nci;
970 :     double
971 :     *omgi = REAL(VECTOR_ELT(Omega, i));
972 :     if (unc) {
973 :     double
974 :     *tmp = Calloc(ncisq, double),
975 :     diagj, one = 1., zero = 0.;
976 :    
977 : bates 429 AZERO(omgi, ncisq);
978 : bates 255 for (j = 0; j < nci; j++) {
979 :     tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
980 :     for (k = j + 1; k < nci; k++) {
981 :     tmp[k*nci + j] = cc[odind++] * diagj;
982 :     }
983 :     }
984 :     F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
985 :     tmp, &nci, &zero, omgi, &nci);
986 :     Free(tmp);
987 :     } else {
988 :     for (j = 0; j < nci; j++) {
989 :     omgi[j * ncip1] = cc[cind++];
990 :     for (k = j + 1; k < nci; k++) {
991 :     omgi[k*nci + j] = cc[odind++];
992 :     }
993 :     }
994 :     }
995 :     cind = odind;
996 :     }
997 :     }
998 :     status[0] = status[1] = 0;
999 :     return x;
1000 :     }
1001 :    
1002 :     /**
1003 :     * Extract the conditional estimates of the fixed effects
1004 :     *
1005 :     * @param x Pointer to an lme object
1006 :     *
1007 :     * @return a numeric vector containing the conditional estimates of
1008 :     * the fixed effects
1009 :     */
1010 : bates 415 SEXP lmer_fixef(SEXP x)
1011 : bates 255 {
1012 :     SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym),
1013 :     cnames = GET_SLOT(x, Matrix_cnamesSym);
1014 : bates 415 int j, pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1];
1015 : bates 255 SEXP val = PROTECT(allocVector(REALSXP, pp1));
1016 :     double
1017 :     *beta = REAL(val),
1018 :     nryyinv; /* negative ryy-inverse */
1019 :    
1020 : bates 415 lmer_invert(x);
1021 : bates 255 Memcpy(beta, REAL(RXXsl) + pp1 * (pp1 - 1), pp1);
1022 :     nryyinv = -REAL(RXXsl)[pp1*pp1 - 1];
1023 :     for (j = 0; j < pp1; j++) beta[j] /= nryyinv;
1024 : bates 415 setAttrib(val, R_NamesSymbol,
1025 :     duplicate(VECTOR_ELT(cnames, length(cnames) - 1)));
1026 : bates 255 UNPROTECT(1);
1027 :     return val;
1028 :     }
1029 :    
1030 :     /**
1031 :     * Extract the conditional modes of the random effects.
1032 :     *
1033 :     * @param x Pointer to an lme object
1034 :     *
1035 :     * @return a vector containing the conditional modes of the random effects
1036 :     */
1037 : bates 415 SEXP lmer_ranef(SEXP x)
1038 : bates 255 {
1039 : bates 429 SEXP RZXP = GET_SLOT(x, Matrix_RZXSym),
1040 : bates 255 cnames = GET_SLOT(x, Matrix_cnamesSym),
1041 : bates 415 flist = GET_SLOT(x, Matrix_flistSym);
1042 :     int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1043 : bates 429 *dims = INTEGER(getAttrib(RZXP, R_DimSymbol)),
1044 : bates 255 *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1045 :     i, ii, jj,
1046 : bates 415 nf = length(flist);
1047 : bates 255 SEXP val = PROTECT(allocVector(VECSXP, nf));
1048 :     double
1049 : bates 429 *b = REAL(RZXP) + dims[0] * (dims[1] - 1),
1050 : bates 255 nryyinv; /* negative ryy-inverse */
1051 :    
1052 : bates 415 lmer_invert(x);
1053 : bates 255 setAttrib(val, R_NamesSymbol,
1054 : bates 415 duplicate(getAttrib(flist, R_NamesSymbol)));
1055 : bates 255 nryyinv = -REAL(GET_SLOT(x, Matrix_RXXSym))[dims[1] * dims[1] - 1];
1056 :     for (i = 0; i < nf; i++) {
1057 : bates 415 SEXP nms, rnms = getAttrib(VECTOR_ELT(flist, i), R_LevelsSymbol);
1058 :     int nci = nc[i], mi = length(rnms);
1059 :     double *bi = b + Gp[i], *mm;
1060 : bates 255
1061 :     SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, mi, nci));
1062 :     setAttrib(VECTOR_ELT(val, i), R_DimNamesSymbol, allocVector(VECSXP, 2));
1063 :     nms = getAttrib(VECTOR_ELT(val, i), R_DimNamesSymbol);
1064 :     SET_VECTOR_ELT(nms, 0, duplicate(rnms));
1065 :     SET_VECTOR_ELT(nms, 1, duplicate(VECTOR_ELT(cnames, i)));
1066 :     mm = REAL(VECTOR_ELT(val, i));
1067 :     for (jj = 0; jj < nci; jj++)
1068 :     for(ii = 0; ii < mi; ii++)
1069 : bates 415 mm[ii + jj * mi] = bi[jj + ii * nci]/nryyinv;
1070 : bates 255 }
1071 :     UNPROTECT(1);
1072 :     return val;
1073 :     }
1074 :    
1075 :     /**
1076 :     * Fill in four symmetric matrices for each level, providing the
1077 :     * information to generate the gradient or the ECME step. The four
1078 :     * matrices are
1079 :     * 1) -m_i\bOmega_i^{-1}
1080 :     * 2) \bB_i\bB_i\trans
1081 :     * 3) \tr\left[\der_{\bOmega_i}\bOmega\left(\bZ\trans\bZ+\bOmega\right)\inv\right]
1082 :     * 4) The term added to 3) to get \tr\left[\der_{\bOmega_i}\bOmega\vb\right]
1083 :     *
1084 :     * @param x pointer to an lme object
1085 :     * @param val pointer to a list of matrices of the correct sizes
1086 :     *
1087 :     * @return val
1088 :     */
1089 :     static
1090 : bates 415 SEXP lmer_firstDer(SEXP x, SEXP val)
1091 : bates 255 {
1092 : bates 429 SEXP bVarP = GET_SLOT(x, Matrix_bVarSym),
1093 :     OmegaP = GET_SLOT(x, Matrix_OmegaSym),
1094 :     RZXP = GET_SLOT(x, Matrix_RZXSym);
1095 :     int *dims = INTEGER(getAttrib(RZXP, R_DimSymbol)),
1096 : bates 415 *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1097 : bates 429 i, nf = length(OmegaP), p = dims[1] - 1;
1098 :     double *RZX = REAL(RZXP),
1099 :     *b = REAL(RZXP) + dims[0] * p;
1100 : bates 255
1101 : bates 415 lmer_invert(x);
1102 : bates 437 /* FIXME: Why is this loop run backwards? It appears it could run forwards. */
1103 : bates 425 for (i = nf - 1; i >= 0; i--) {
1104 : bates 429 SEXP bVPi = VECTOR_ELT(bVarP, i);
1105 :     int *ddims = INTEGER(getAttrib(bVPi, R_DimSymbol)), j, k;
1106 :     int nci = ddims[0];
1107 : bates 415 int ncisqr = nci * nci, RZXrows = Gp[i + 1] - Gp[i];
1108 : bates 429 int nlev = RZXrows/nci;
1109 :     double *RZXi = RZX + Gp[i], *bVi = REAL(bVPi),
1110 : bates 422 *bi = b + Gp[i], *mm = REAL(VECTOR_ELT(val, i)),
1111 :     *tmp = Memcpy(Calloc(ncisqr, double),
1112 : bates 429 REAL(VECTOR_ELT(OmegaP, i)), ncisqr),
1113 : bates 422 dlev = (double) nlev,
1114 :     one = 1., zero = 0.;
1115 :    
1116 : bates 255 if (nci == 1) {
1117 :     int ione = 1;
1118 : bates 429 mm[0] = ((double) nlev)/tmp[0];
1119 : bates 415 mm[1] = F77_CALL(ddot)(&nlev, bi, &ione, bi, &ione);
1120 : bates 429 mm[2] = 0.;
1121 :     for (k = 0; k < nlev; k++) mm[2] += bVi[k];
1122 : bates 424 mm[3] = 0.;
1123 : bates 255 for (j = 0; j < p; j++) {
1124 : bates 415 mm[3] += F77_CALL(ddot)(&RZXrows, RZXi + j * dims[0], &ione,
1125 :     RZXi + j * dims[0], &ione);
1126 : bates 255 }
1127 :     } else {
1128 : bates 429 AZERO(mm, 4 * ncisqr);
1129 : bates 255 F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1130 :     if (j)
1131 :     error("Omega[[%d]] is not positive definite", i + 1);
1132 :     F77_CALL(dtrtri)("U", "N", &nci, tmp, &nci, &j);
1133 :     if (j)
1134 :     error("Omega[[%d]] is not positive definite", i + 1);
1135 :     F77_CALL(dsyrk)("U", "N", &nci, &nci, &dlev, tmp, &nci,
1136 :     &zero, mm, &nci);
1137 :     mm += ncisqr; /* \bB_i term */
1138 : bates 415 F77_CALL(dsyrk)("U", "N", &nci, &nlev, &one, bi, &nci,
1139 : bates 255 &zero, mm, &nci);
1140 : bates 429 mm += ncisqr; /* Sum of diagonal blocks of the inverse
1141 :     * (Z'Z+Omega)^{-1} */
1142 :     for (j = 0; j < ncisqr; j++) {
1143 :     for (k = 0; k < nlev; k++) mm[j] += bVi[j + k*ncisqr];
1144 :     }
1145 : bates 255 mm += ncisqr; /* Extra term for \vb */
1146 :     for (j = 0; j < p; j++) {
1147 :     F77_CALL(dsyrk)("U", "N", &nci, &nlev, &one,
1148 : bates 415 RZXi + j * dims[0], &nci,
1149 : bates 255 &one, mm, &nci);
1150 :     }
1151 :     }
1152 : bates 424 Free(tmp);
1153 : bates 255 }
1154 :     return val;
1155 :     }
1156 :    
1157 :     /**
1158 : bates 422 * Return a length nf list of arrays of dimension (nci, nci, 4). The
1159 :     * values of these arrays are assigned in lmer_firstDer.
1160 : bates 255 *
1161 :     * @param nf number of factors
1162 :     * @param nc vector of number of columns per factor
1163 :     *
1164 :     * @return pointer to a list of REAL arrays
1165 :     */
1166 :     static
1167 :     SEXP EM_grad_array(int nf, const int nc[])
1168 :     {
1169 : bates 424 SEXP val = PROTECT(allocVector(VECSXP, nf));
1170 : bates 422 int i;
1171 : bates 255
1172 :     for (i = 0; i < nf; i++) {
1173 : bates 422 SET_VECTOR_ELT(val, i, alloc3Darray(REALSXP, nc[i], nc[i], 4));
1174 : bates 255 }
1175 : bates 422 UNPROTECT(1);
1176 : bates 255 return val;
1177 :     }
1178 :    
1179 :     /**
1180 :     * Fill in the 4-dimensional vector of linear combinations of the
1181 : bates 415 * firstDer array according to whether ECME steps or the gradient are
1182 : bates 255 * needed and to whether or not REML is being used.
1183 :     *
1184 :     * @param cc coefficient vector to be filled in
1185 :     * @param EM non-zero for ECME steps, zero for gradient
1186 :     * @param REML non-zero for REML, zero for ML
1187 :     * @param ns ns[0] is p+1, ns[1] is n
1188 :     *
1189 :     * @return cc with the coefficients filled in
1190 :     */
1191 :     static
1192 :     double *EM_grad_lc(double *cc, int EM, int REML, int ns[])
1193 :     {
1194 :     cc[0] = EM ? 0. : -1.;
1195 :     cc[1] = (double)(ns[1] - (REML ? ns[0] - 1 : 0));
1196 :     cc[2] = 1.;
1197 :     cc[3] = REML ? 1. : 0.;
1198 :     return cc;
1199 :     }
1200 :    
1201 :    
1202 :     /**
1203 :     * Print the verbose output in the ECME iterations
1204 :     *
1205 :     * @param x pointer to an ssclme object
1206 :     * @param iter iteration number
1207 :     * @param REML non-zero for REML, zero for ML
1208 : bates 429 * @param firstDer arrays for calculating ECME steps and the first derivative
1209 :     * @param val Pointer to a list of arrays to receive the calculated values
1210 : bates 255 */
1211 :     static
1212 :     void EMsteps_verbose_print(SEXP x, int iter, int REML, SEXP firstDer, SEXP val)
1213 :     {
1214 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym),
1215 :     pMat = VECTOR_ELT(val, 2);
1216 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1217 :     *Its = INTEGER(VECTOR_ELT(val, 0)),
1218 :     i, ifour = 4, ii, ione = 1, jj, nf = length(Omega),
1219 :     niter = INTEGER(getAttrib(pMat, R_DimSymbol))[0];
1220 :     double
1221 :     *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),
1222 :     *cc = EM_grad_lc(Calloc(4, double), 0, REML, nc + nf),
1223 :     *Devs = REAL(VECTOR_ELT(val, 1)),
1224 :     *pars = REAL(pMat) + iter,
1225 :     *grds = REAL(VECTOR_ELT(val, 3)) + iter,
1226 :     one = 1., zero = 0.;
1227 :    
1228 : bates 415 lmer_factor(x);
1229 : bates 255 if (iter == 0) Rprintf(" EM iterations\n");
1230 :     Rprintf("%3d %.3f", Its[iter] = iter, Devs[iter] = dev[REML ? 1 : 0]);
1231 :     for (i = 0; i < nf; i++) {
1232 :     int nci = nc[i], ncip1 = nci + 1, ncisqr = nci * nci;
1233 :     double
1234 :     *Omgi = REAL(VECTOR_ELT(Omega, i)),
1235 :     *Grad = Calloc(ncisqr, double);
1236 :    
1237 :     /* diagonals */
1238 : bates 443 Rprintf(" (%#8g", *pars = Omgi[0]);
1239 :     pars += niter;
1240 :     for (jj = 1; jj < nci; jj++, pars += niter) {
1241 : bates 255 Rprintf(" %#8g", *pars = Omgi[jj * ncip1]);
1242 :     }
1243 :     for (jj = 1; jj < nci; jj++) /* offdiagonals */
1244 :     for (ii = 0; ii < jj; ii++, pars += niter)
1245 :     Rprintf(" %#8g", *pars = Omgi[ii + jj * nci]);
1246 :     /* Evaluate and print the gradient */
1247 :     F77_CALL(dgemv)("N", &ncisqr, &ifour, &one,
1248 :     REAL(VECTOR_ELT(firstDer, i)), &ncisqr,
1249 :     cc, &ione, &zero, Grad, &ione);
1250 : bates 443 Rprintf(":%#8.3g", *grds = Grad[0]);
1251 :     grds += niter;
1252 : bates 255 /* diagonals */
1253 : bates 443 for (jj = 1; jj < nci; jj++, grds += niter) {
1254 :     Rprintf(" %#8.3g", *grds = Grad[jj * ncip1]);
1255 : bates 255 }
1256 :     for (jj = 1; jj < nci; jj++) /* offdiagonals */
1257 :     for (ii = 0; ii < jj; ii++, grds += niter)
1258 : bates 443 Rprintf(" %#8.3g", *grds = Grad[ii + jj * nci]);
1259 :     Rprintf(")");
1260 : bates 255 Free(Grad);
1261 :     }
1262 :     Rprintf("\n");
1263 :     Free(cc);
1264 :     }
1265 :    
1266 :     /**
1267 :     * Perform ECME steps for the REML or ML criterion.
1268 :     *
1269 :     * @param x pointer to an ssclme object
1270 :     * @param nsteps pointer to an integer scalar - the number of ECME steps to perform
1271 :     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1272 :     * @param verb pointer to a logical scalar indicating verbose output
1273 :     *
1274 :     * @return R_NilValue if verb == FALSE, otherwise a list of iteration
1275 :     *numbers, deviances, parameters, and gradients.
1276 :     */
1277 : bates 415 SEXP lmer_ECMEsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP Verbp)
1278 : bates 255 {
1279 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym),
1280 : bates 415 flist = GET_SLOT(x, Matrix_flistSym),
1281 : bates 255 val = R_NilValue;
1282 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1283 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1284 :     REML = asLogical(REMLp),
1285 :     i, ifour = 4, info, ione = 1, iter,
1286 :     nEM = asInteger(nsteps),
1287 :     nf = length(Omega),
1288 :     verb = asLogical(Verbp);
1289 :     double
1290 :     *cc = EM_grad_lc(Calloc(4, double), 1, REML, nc + nf),
1291 :     zero = 0.0;
1292 :     SEXP firstDer = PROTECT(EM_grad_array(nf, nc));
1293 :    
1294 : bates 429 lmer_firstDer(x, firstDer);
1295 : bates 255 if (verb) {
1296 :     int nEMp1 = nEM + 1, npar = coef_length(nf, nc);
1297 :     val = PROTECT(allocVector(VECSXP, 4));
1298 :     SET_VECTOR_ELT(val, 0, allocVector(INTSXP, nEMp1));
1299 :     SET_VECTOR_ELT(val, 1, allocVector(REALSXP, nEMp1));
1300 :     SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, nEMp1, npar));
1301 :     SET_VECTOR_ELT(val, 3, allocMatrix(REALSXP, nEMp1, npar));
1302 : bates 429 EMsteps_verbose_print(x, 0, REML, firstDer, val);
1303 : bates 255 }
1304 :     for (iter = 0; iter < nEM; iter++) {
1305 :     for (i = 0; i < nf; i++) {
1306 :     int nci = nc[i], ncisqr = nci * nci;
1307 :     double *Omgi = REAL(VECTOR_ELT(Omega, i)),
1308 : bates 415 mult = 1./
1309 :     ((double) length(getAttrib(VECTOR_ELT(flist, i),
1310 :     R_LevelsSymbol)));
1311 : bates 255
1312 :     F77_CALL(dgemm)("N", "N", &ncisqr, &ione, &ifour, &mult,
1313 :     REAL(VECTOR_ELT(firstDer, i)), &ncisqr,
1314 :     cc, &ifour, &zero, Omgi, &ncisqr);
1315 :     F77_CALL(dpotrf)("U", &nci, Omgi, &nci, &info);
1316 :     if (info)
1317 :     error("DPOTRF in ECME update gave code %d", info);
1318 :     F77_CALL(dpotri)("U", &nci, Omgi, &nci, &info);
1319 :     if (info)
1320 :     error("Matrix inverse in ECME update gave code %d", info);
1321 :     }
1322 :     status[0] = status[1] = 0;
1323 : bates 429 lmer_firstDer(x, firstDer);
1324 : bates 255 if (verb) EMsteps_verbose_print(x, iter + 1, REML, firstDer, val);
1325 :     }
1326 : bates 415 lmer_factor(x);
1327 : bates 255 if (verb) UNPROTECT(1);
1328 :     UNPROTECT(1);
1329 :     return val;
1330 :     }
1331 :    
1332 : bates 415 SEXP lmer_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1333 : bates 255 {
1334 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1335 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1336 :     dind, i, ifour = 4, info, ione = 1, nf = length(Omega),
1337 :     odind, unc = asLogical(Uncp);
1338 :     SEXP
1339 : bates 415 firstDer = lmer_firstDer(x, PROTECT(EM_grad_array(nf, nc))),
1340 : bates 255 val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1341 :     double
1342 :     *cc = EM_grad_lc(Calloc(4, double), 0,
1343 :     asInteger(REMLp), nc + nf),
1344 :     one = 1.0, zero = 0.0;
1345 :    
1346 :     dind = 0; /* index into val for diagonals */
1347 :     for (i = 0; i < nf; i++) {
1348 :     int nci = nc[i], ncisqr = nci * nci;
1349 :     double
1350 :     *Omgi = REAL(VECTOR_ELT(Omega, i)),
1351 :     *tmp = Calloc(ncisqr, double);
1352 :    
1353 :     F77_CALL(dgemm)("N", "N", &ncisqr, &ione, &ifour, &one,
1354 :     REAL(VECTOR_ELT(firstDer, i)), &ncisqr,
1355 :     cc, &ifour, &zero, tmp, &ncisqr);
1356 :     if (nci == 1) {
1357 :     REAL(val)[dind++] = (unc ? Omgi[0] : 1.) * tmp[0];
1358 :     } else {
1359 :     int ii, j, ncip1 = nci + 1;
1360 :    
1361 :     odind = dind + nci; /* index into val for off-diagonals */
1362 :     if (unc) {
1363 :     double *chol = Memcpy(Calloc(ncisqr, double),
1364 :     REAL(VECTOR_ELT(Omega, i)), ncisqr),
1365 :     *tmp2 = Calloc(ncisqr, double);
1366 :    
1367 :     /* Overwrite the gradient with respect to positions in
1368 :     * Omega[[i]] by the gradient with respect to the
1369 :     * unconstrained parameters.*/
1370 :    
1371 :     F77_CALL(dpotrf)("U", &nci, chol, &nci, &info);
1372 :     if (info)
1373 :     error("Omega[[%d]] is not positive definite", i + 1);
1374 :     /* tmp2 := chol %*% tmp using only upper triangle of tmp */
1375 :     F77_CALL(dsymm)("R", "U", &nci, &nci, &one, tmp, &nci,
1376 :     chol, &nci, &zero, tmp2, &nci);
1377 :     /* full symmetric product gives diagonals */
1378 :     F77_CALL(dtrmm)("R", "U", "T", "N", &nci, &nci, &one, chol, &nci,
1379 :     Memcpy(tmp, tmp2, ncisqr), &nci);
1380 :     /* overwrite upper triangle with gradients for positions in L' */
1381 :     for (ii = 1; ii < nci; ii++) {
1382 :     for (j = 0; j < ii; j++) {
1383 :     tmp[j + ii*nci] = chol[j*ncip1] * tmp2[j + ii*nci];
1384 :     tmp[ii + j*nci] = 0.;
1385 :     }
1386 :     }
1387 :     Free(chol); Free(tmp2);
1388 :     }
1389 :     for (j = 0; j < nci; j++) {
1390 :     REAL(val)[dind + j] = tmp[j * ncip1];
1391 :     for (ii = 0; ii < j; ii++) /* offdiagonals count twice */
1392 :     REAL(val)[odind++] = 2. * tmp[ii + j * nci];
1393 :     }
1394 :     dind = odind;
1395 :     }
1396 :     Free(tmp);
1397 :     }
1398 :     UNPROTECT(2);
1399 :     Free(cc);
1400 :     return val;
1401 :     }
1402 :    
1403 :     /**
1404 :     * Fill in five symmetric matrices, providing the
1405 :     * information to generate the Hessian.
1406 :    
1407 :     * @param x pointer to an lme object
1408 :     * @param val ignored at present
1409 :     *
1410 :     * @return val an array consisting of five symmetric faces
1411 :     */
1412 : bates 296 static
1413 : bates 415 SEXP lmer_secondDer(SEXP x, SEXP Valp)
1414 : bates 255 {
1415 :     SEXP
1416 :     D = GET_SLOT(x, Matrix_DSym),
1417 :     Omega = GET_SLOT(x, Matrix_OmegaSym),
1418 : bates 429 RZXP = GET_SLOT(x, Matrix_RZXSym),
1419 : bates 255 levels = GET_SLOT(x, R_LevelsSymbol),
1420 :     val;
1421 : bates 429 int *dRZX = INTEGER(getAttrib(RZXP, R_DimSymbol)),
1422 : bates 255 *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1423 :     Q, Qsqr, RZXpos, facepos,
1424 :     i, ione = 1, j, nf = length(Omega), p = dRZX[1] - 1, pos;
1425 :     SEXP
1426 : bates 415 firstDer = lmer_firstDer(x, PROTECT(EM_grad_array(nf, nc)));
1427 : bates 255 double
1428 : bates 429 *RZX = REAL(RZXP),
1429 :     *b = REAL(RZXP) + dRZX[0] * p,
1430 : bates 255 *bbface, /* vec of second faces of firstDer elts */
1431 :     one = 1.,
1432 :     zero = 0.;
1433 :    
1434 :     Q = 0; /* number of rows and columns in the result */
1435 :     for (i = 0; i < nf; i++) Q += nc[i] * nc[i];
1436 :     Qsqr = Q * Q;
1437 :     bbface = Calloc(Q, double);
1438 : bates 275 val = PROTECT(alloc3Darray(REALSXP, Q, Q, 5));
1439 : bates 429 AZERO(REAL(val), Qsqr * 5);
1440 : bates 255
1441 :     pos = 0;
1442 :     for (i = 0; i < nf; i++) {
1443 :     int nci = nc[i], ncisqr = nci * nci;
1444 :     double *fDi = REAL(VECTOR_ELT(firstDer, i)),
1445 :     mult = 1./((double) length(VECTOR_ELT(levels, i)));
1446 :    
1447 :     Memcpy(bbface + pos, fDi + ncisqr, ncisqr);
1448 :     /* outer product of the third face of firstDer on the diagonal
1449 :     * of the third face of val */
1450 :     F77_CALL(dsyr)("U", &ncisqr, &mult, fDi + 2 * ncisqr, &ione,
1451 :     REAL(val) + 2 * Qsqr + pos * Q, &Q);
1452 :     pos += ncisqr;
1453 :     }
1454 :     /* fifth face of val is outer product of bbface */
1455 :     F77_CALL(dsyr)("U", &Q, &one, bbface, &ione, REAL(val) + 4 * Qsqr, &Q);
1456 :     /* fourth face from \bb\trans\der\vb\der\bb */
1457 : bates 429 AZERO(REAL(val) + 3 * Qsqr, Qsqr); /* zero accumulator */
1458 : bates 255 RZXpos = 0;
1459 :     facepos = 0;
1460 :     for (i = 0; i < nf; i++) {
1461 :     int ii, jj, nci = nc[i], ncisqr = nci * nci, nctp = nci * p,
1462 :     nlev = length(VECTOR_ELT(levels, i));
1463 :     int maxpq = (p > nci) ? p : nci;
1464 :     double
1465 :     *Di = REAL(VECTOR_ELT(D, i)),
1466 :     *arr = Calloc(ncisqr * maxpq, double), /* tmp 3Darray */
1467 :     *face = REAL(val) + 3 * Qsqr,
1468 :     *mat = Calloc(nci * maxpq, double); /* tmp matrix */
1469 :    
1470 :     for (j = 0; j < nlev; j++) {
1471 :     F77_CALL(dgemm)("T", "T", &p, &nci, &nci,
1472 :     &one, RZX + j * nci, dRZX, Di + j * ncisqr, &nci,
1473 :     &zero, mat, &p);
1474 :     F77_CALL(dgemm)("N", "N", &nctp, &nci, &ione,
1475 :     &one, mat, &nctp, b + j * nci, &ione,
1476 :     &zero, arr, &nctp);
1477 :     F77_CALL(dsyrk)("U", "T", &ncisqr, &p, &one, arr, &p,
1478 :     &one, face + facepos, &Q);
1479 :     /* Add the D_{i,j}^{-T/2} term */
1480 :     Memcpy(mat, Di + j * ncisqr, ncisqr);
1481 :     for (jj = 1; jj < nci; jj++) { /* transpose mat */
1482 :     for (ii = 0; ii < jj; ii++) {
1483 :     mat[jj + ii * nci] = mat[ii + jj * nci];
1484 :     mat[ii + jj * nci] = 0.;
1485 :     }
1486 :     }
1487 :     F77_CALL(dgemm)("N", "N", &ncisqr, &nci, &ione,
1488 :     &one, mat, &ncisqr, b + j * nci, &ione,
1489 :     &zero, arr, &ncisqr);
1490 :     /* FIXME: Next call could be dsyr (it's rank one). */
1491 :     F77_CALL(dsyrk)("U", "T", &ncisqr, &nci, &one, arr, &nci,
1492 :     &one, face + facepos, &Q);
1493 :    
1494 :     }
1495 :     RZXpos += nci * nlev;
1496 :     facepos += ncisqr;
1497 :     Free(arr); Free(mat);
1498 :     }
1499 :     UNPROTECT(2);
1500 :     Free(bbface);
1501 :     return val;
1502 :     }
1503 :    
1504 : bates 274 /**
1505 :     * Return the unscaled variances
1506 :     *
1507 : bates 415 * @param x pointer to an lmer object
1508 : bates 274 *
1509 :     * @return a list similar to the Omega list with the unscaled variances
1510 :     */
1511 : bates 415 SEXP lmer_variances(SEXP x)
1512 : bates 274 {
1513 :     SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1514 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1515 :     i, nf = length(Omg);
1516 :    
1517 :     for (i = 0; i < nf; i++) {
1518 :     double *mm = REAL(VECTOR_ELT(Omg, i));
1519 :     int j, nci = nc[i];
1520 :    
1521 :     F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1522 :     if (j) /* shouldn't happen */
1523 :     error("DPOTRF returned error code %d on Omega[%d]",
1524 :     j, i + 1);
1525 :     F77_CALL(dpotri)("U", &nci, mm, &nci, &j);
1526 :     if (j) /* shouldn't happen */
1527 :     error("DTRTRI returned error code %d on Omega[%d]",
1528 :     j, i + 1);
1529 :     nlme_symmetrize(mm, nci);
1530 :     }
1531 :     UNPROTECT(1);
1532 :     return Omg;
1533 :     }
1534 : bates 410
1535 : bates 466 SEXP lmer_Crosstab(SEXP flist)
1536 :     {
1537 :     SEXP val;
1538 :     int i, nf = length(flist), nobs;
1539 :     int *nc = Calloc(nf, int);
1540 :    
1541 :     if (!(nf > 0 && isNewList(flist)))
1542 :     error("flist must be a non-empty list");
1543 :     nobs = length(VECTOR_ELT(flist, 0));
1544 :     if (nobs < 1) error("flist[[0]] must be a non-null factor");
1545 :     for (i = 0; i < nf; i++) {
1546 :     SEXP fi = VECTOR_ELT(flist, i);
1547 :     if (!(isFactor(fi) && length(fi) == nobs))
1548 :     error("flist[[%d]] must be a factor of length %d",
1549 :     i + 1, nobs);
1550 :     nc[i] = 1;
1551 :     }
1552 :     val = lmer_crosstab(flist, nobs, nc);
1553 :     Free(nc);
1554 :     return val;
1555 :     }

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