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

1 : bates 10 #include "ssclme.h"
2 :    
3 : bates 176 #define slot_dup(dest, src, sym) SET_SLOT(dest, sym, duplicate(GET_SLOT(src, sym)))
4 :    
5 :     /**
6 :     * Using the sscCrosstab object from the grouping factors, generate
7 :     * the slots in an ssclme object related to the symmetric sparse
8 :     * matrix representation of Z'Z. If the model matrices for the
9 :     * grouping factors have only one column each then the structure can
10 :     * be copied, otherwise it must be generated from the sscCrosstab and
11 :     * the number of columns per grouping factor.
12 :     *
13 :     * @param nf number of factors
14 :     * @param nc vector of length nf+2 with number of columns in model matrices
15 :     * @param ctab pointer to the sscCrosstab object
16 :     * @param ssc pointer to an ssclme object to be filled out
17 :     */
18 : bates 10 static
19 :     void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)
20 :     {
21 :     int *snc, i, copyonly = 1;
22 :    
23 : bates 176 SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));
24 :     snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
25 :     for (i = 0; i <= nf; i++) {
26 :     snc[i] = nc[i];
27 :     if (nc[i] > 1 && i < nf) copyonly = 0;
28 : bates 10 }
29 :     if (copyonly) {
30 : bates 176 slot_dup(ssc, ctab, Matrix_pSym);
31 :     slot_dup(ssc, ctab, Matrix_iSym);
32 :     slot_dup(ssc, ctab, Matrix_xSym);
33 :     slot_dup(ssc, ctab, Matrix_DimSym);
34 :     slot_dup(ssc, ctab, Matrix_GpSym);
35 :     return;
36 :     }
37 :     {
38 : bates 10 int
39 :     *GpIn = INTEGER(GET_SLOT(ctab, Matrix_GpSym)),
40 :     *GpOut,
41 :     *AiIn = INTEGER(GET_SLOT(ctab, Matrix_iSym)),
42 :     *AiOut,
43 :     *ApIn = INTEGER(GET_SLOT(ctab, Matrix_pSym)),
44 :     *ApOut,
45 :     nIn = GpIn[nf], nOut, nzOut,
46 :     *dims,
47 :     *map = Calloc(nIn + 1, int), /* column number map */
48 :     *ncc = Calloc(nIn, int); /* number of columns out for this
49 :     * col in */
50 :    
51 :     SET_SLOT(ssc, Matrix_GpSym, allocVector(INTSXP, nf + 1));
52 :     GpOut = INTEGER(GET_SLOT(ssc, Matrix_GpSym));
53 :     map[0] = GpOut[0] = 0;
54 :     for (i = 0; i < nf; i++) {
55 :     int j, GpIni = GpIn[i], GpInip1 = GpIn[i+1], nci = nc[i];
56 :     GpOut[i+1] = GpOut[i] + (GpInip1 - GpIni) * nci;
57 :     for (j = GpIni; j < GpInip1; j++) {
58 :     ncc[j] = nci;
59 :     map[j+1] = map[j] + nci;
60 :     }
61 :     }
62 :     nOut = GpOut[nf]; /* size of output matrix */
63 : bates 176 SET_SLOT(ssc, Matrix_DimSym, allocVector(INTSXP, 2));
64 : bates 10 dims = INTEGER(GET_SLOT(ssc, Matrix_DimSym));
65 :     dims[0] = dims[1] = nOut;
66 :     SET_SLOT(ssc, Matrix_pSym, allocVector(INTSXP, nOut + 1));
67 :     ApOut = INTEGER(GET_SLOT(ssc, Matrix_pSym));
68 :     ApOut[0] = 0;
69 :     for (i = 0; i < nf; i++) { /* determine the column pointers */
70 :     int j, jout = GpOut[i], nci = nc[i], p2 = GpIn[i+1];
71 :     for (j = GpIn[i]; j < p2; j++) {
72 :     int k, nj = 0, p3 = ApIn[j+1];
73 :     for (k = ApIn[j]; k < p3; k++) {
74 :     nj += ncc[AiIn[k]];
75 :     }
76 :     nj -= nci - 1;
77 :     ApOut[jout+1] = ApOut[jout] + nj;
78 :     jout++;
79 :     for (k = 1; k < nci; k++) {
80 :     ApOut[jout+1] = ApOut[jout] + nj + k;
81 :     jout++;
82 :     }
83 :     }
84 :     }
85 :     nzOut = ApOut[nOut]; /* number of non-zeros in output */
86 :     SET_SLOT(ssc, Matrix_xSym, allocVector(REALSXP, nzOut));
87 :     memset(REAL(GET_SLOT(ssc, Matrix_xSym)), 0,
88 :     sizeof(double) * nzOut);
89 :     SET_SLOT(ssc, Matrix_iSym, allocVector(INTSXP, nzOut));
90 :     AiOut = INTEGER(GET_SLOT(ssc, Matrix_iSym));
91 :     for (i = 0; i < nf; i++) { /* fill in the rows */
92 :     int j, jj, nci = nc[i], p2 = GpIn[i+1];
93 :     for (j = GpIn[i]; j < p2; j++) { /* first col for input col */
94 : bates 168 int k, mj = map[j], p3 = ApIn[j+1], pos = ApOut[mj];
95 :     for (k = ApIn[j]; k < p3; k++) {
96 :     int ii = AiIn[k], ncci = ncc[ii];
97 :     AiOut[pos++] = map[ii];
98 :     if (ii < j) { /* above the diagonal */
99 :     for (jj = 1; jj < ncci; jj++) {
100 :     AiOut[pos+1] = AiOut[pos] + 1;
101 :     pos++;
102 :     }
103 : bates 10 }
104 : bates 168 for (jj = 1; jj < nci; jj++) { /* repeat the column adding
105 :     * another diagonal element */
106 :     int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];
107 :     Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);
108 : bates 176 AiOut[ApOut[mjj + 1] - 1] = mjj;
109 : bates 168 }
110 : bates 10 }
111 :     }
112 :     }
113 :     Free(map); Free(ncc);
114 :     }
115 :     }
116 :    
117 : bates 159 /**
118 :     * Calculate and store the maximum number of off-diagonal elements in
119 :     * the inverse of L, based on the elimination tree. The maximum is
120 :     * itself stored in the Parent array. (FIXME: come up with a better design.)
121 :     *
122 :     * @param n number of columns in the matrix
123 :     * @param Parent elimination tree for the matrix
124 :     */
125 :     static void ssclme_calc_maxod(int n, int Parent[])
126 : bates 10 {
127 : bates 159 int *sz = Calloc(n, int), i, mm = -1;
128 : bates 10 for (i = n - 1; i >= 0; i--) {
129 : bates 159 sz[i] = (Parent[i] < 0) ? 0 : (1 + sz[Parent[i]]);
130 :     if (sz[i] > mm) mm = sz[i];
131 : bates 10 }
132 : bates 159 Parent[n] = mm;
133 : bates 10 Free(sz);
134 :     }
135 :    
136 : bates 176 /**
137 :     * Create an ssclme object from a list of grouping factors, sorted in
138 :     * order of non-increasing numbers of levels, and an integer vector of
139 :     * the number of columns in the model matrices. There is one more
140 :     * element in ncv than in facs. The last element is the number of
141 :     * columns in the model matrix for the fixed effects plus the
142 :     * response. (i.e. p+1)
143 :     *
144 :     * @param facs pointer to a list of grouping factors
145 :     * @param ncv pointer to an integer vector of number of columns per model matrix
146 :     *
147 :     * @return pointer to an ssclme object
148 :     */
149 : bates 10 SEXP
150 : bates 168 ssclme_create(SEXP facs, SEXP ncv)
151 : bates 10 {
152 :     SEXP ctab, nms, ssc, tmp,
153 : bates 21 val = PROTECT(allocVector(VECSXP, 2)),
154 :     dd = PROTECT(allocVector(INTSXP, 3)); /* dimensions of 3-D arrays */
155 : bates 159 int *Ai, *Ap, *Gp, *Lp, *Parent,
156 : bates 21 *nc, Lnz, i, nf = length(facs), nzcol, pp1,
157 :     *dims = INTEGER(dd);
158 : bates 10
159 :     if (length(ncv) != (nf + 1))
160 :     error("length of nc (%d) should be length of facs (%d) + 1",
161 :     length(ncv), nf);
162 :     SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS("ssclme")));
163 :     ssc = VECTOR_ELT(val, 0);
164 :     /* Pairwise cross-tabulation */
165 :     ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));
166 : bates 143 SET_VECTOR_ELT(val, 1, sscCrosstab_groupedPerm(ctab));
167 : bates 10 if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */
168 :     ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
169 :     1, INTEGER(VECTOR_ELT(val, 1)),
170 :     INTEGER(GET_SLOT(ctab, Matrix_pSym)),
171 :     INTEGER(GET_SLOT(ctab, Matrix_iSym)));
172 :     }
173 :     ssclme_copy_ctab(nf, INTEGER(ncv), ctab, ssc);
174 :     UNPROTECT(1); /* ctab */
175 :    
176 :     nzcol = INTEGER(GET_SLOT(ssc, Matrix_DimSym))[1];
177 :     Gp = INTEGER(GET_SLOT(ssc, Matrix_GpSym));
178 :     Ap = INTEGER(GET_SLOT(ssc, Matrix_pSym));
179 :     Ai = INTEGER(GET_SLOT(ssc, Matrix_iSym));
180 :     nc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
181 :     nc[nf + 1] = length(VECTOR_ELT(facs, 0)); /* number of observations */
182 :     /* Create slots */
183 :     pp1 = nc[nf];
184 :     SET_SLOT(ssc, Matrix_XtXSym, allocMatrix(REALSXP, pp1, pp1));
185 :     SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));
186 :     SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));
187 :     SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));
188 : bates 101 /* Zero symmetric matrices (cosmetic) */
189 : bates 10 memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,
190 :     sizeof(double) * pp1 * pp1);
191 :     memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,
192 :     sizeof(double) * pp1 * pp1);
193 :     SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));
194 :     Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));
195 : bates 159 SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol + 1));
196 : bates 10 Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));
197 :     SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));
198 :     SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));
199 : bates 372 R_ldl_symbolic(nzcol, Ap, Ai, Lp, Parent,
200 : bates 10 (int *) NULL, (int *) NULL); /* P & Pinv */
201 : bates 159 ssclme_calc_maxod(nzcol, Parent);
202 : bates 10 Lnz = Lp[nzcol];
203 :     SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));
204 :     SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));
205 :     SET_SLOT(ssc, Matrix_OmegaSym, allocVector(VECSXP, nf));
206 :     tmp = GET_SLOT(ssc, Matrix_OmegaSym);
207 :     setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));
208 :     for (i = 0; i < nf; i++) {
209 :     SET_VECTOR_ELT(tmp, i, allocMatrix(REALSXP, nc[i], nc[i]));
210 :     memset(REAL(VECTOR_ELT(tmp, i)), 0,
211 :     sizeof(double) * nc[i] * nc[i]);
212 :     }
213 :     SET_SLOT(ssc, Matrix_devianceSym, allocVector(REALSXP, 2));
214 :     tmp = GET_SLOT(ssc, Matrix_devianceSym);
215 :     setAttrib(tmp, R_NamesSymbol, allocVector(STRSXP, 2));
216 :     nms = getAttrib(tmp, R_NamesSymbol);
217 :     SET_STRING_ELT(nms, 0, mkChar("ML"));
218 :     SET_STRING_ELT(nms, 1, mkChar("REML"));
219 :     SET_SLOT(ssc, Matrix_devCompSym, allocVector(REALSXP, 4));
220 :     SET_SLOT(ssc, Matrix_statusSym, allocVector(LGLSXP, 2));
221 :     tmp = GET_SLOT(ssc, Matrix_statusSym);
222 :     LOGICAL(tmp)[0] = LOGICAL(tmp)[1] = 0;
223 :     setAttrib(tmp, R_NamesSymbol, allocVector(STRSXP, 2));
224 :     nms = getAttrib(tmp, R_NamesSymbol);
225 :     SET_STRING_ELT(nms, 0, mkChar("factored"));
226 :     SET_STRING_ELT(nms, 1, mkChar("inverted"));
227 :     SET_SLOT(ssc, Matrix_bVarSym, allocVector(VECSXP, nf));
228 :     tmp = GET_SLOT(ssc, Matrix_bVarSym);
229 :     setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));
230 :     for (i = 0; i < nf; i++) {
231 : bates 21 int nci = nc[i], mi = (Gp[i+1] - Gp[i])/nc[i];
232 : bates 10
233 : bates 21 dims[0] = dims[1] = nci;
234 :     dims[2] = mi;
235 :     SET_VECTOR_ELT(tmp, i, allocArray(REALSXP, dd));
236 : bates 10 memset(REAL(VECTOR_ELT(tmp, i)), 0,
237 : bates 21 sizeof(double) * nci * nci * mi);
238 : bates 10 }
239 : bates 21 UNPROTECT(2);
240 : bates 10 return val;
241 :     }
242 :    
243 : bates 176 /**
244 :     * Copy information on Z'Z accumulated in the bVar array to Z'Z
245 :     *
246 :     * @param ncj number of columns in this block
247 :     * @param Gpj initial column for this group
248 :     * @param Gpjp initial column for the next group
249 :     * @param bVj pointer to the ncj x ncj x mj array to be filled
250 :     * @param Ap column pointer array for Z'Z
251 :     * @param Ai row indices for Z'Z
252 :     * @param Ax elements of Z'Z
253 :     */
254 : bates 10 static
255 :     void bVj_to_A(int ncj, int Gpj, int Gpjp, const double bVj[],
256 :     const int Ap[], const int Ai[], double Ax[])
257 :     {
258 :     int i, diag, k;
259 :     for (i = Gpj; i < Gpjp; i += ncj) {
260 :     for (k = 0; k < ncj; k++) {
261 :     diag = Ap[i + k + 1] - 1;
262 :     if (Ai[diag] != i+k)
263 :     error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",
264 :     diag, i+k, Ai[diag]);
265 :     Memcpy(Ax + diag - k, bVj + (i+k-Gpj)*ncj, k + 1);
266 :     }
267 :     }
268 :     }
269 :    
270 : bates 176 /**
271 :     * Copy the dimnames from the list of grouping factors and the model
272 :     * matrices for the grouping factors into the appropriate parts of the
273 :     * ssclme object.
274 :     *
275 :     * @param x pointer to an ssclme object
276 :     * @param facs pointer to a list of factors
277 :     * @param mmats pointer to a list of model matrices
278 :     *
279 :     * @return NULL
280 :     */
281 : bates 10 SEXP
282 : bates 22 ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)
283 :     {
284 :     SEXP bVar = GET_SLOT(x, Matrix_bVarSym),
285 :     nms2 = PROTECT(allocVector(VECSXP, 2)),
286 :     nms3 = PROTECT(allocVector(VECSXP, 3));
287 :     int i, nf = length(mmats) - 1;
288 :     SEXP xcols = VECTOR_ELT(GetArrayDimnames(VECTOR_ELT(mmats, nf)), 1);
289 :    
290 :     for (i = 0; i < nf; i++) {
291 :     SEXP cnms = VECTOR_ELT(GetArrayDimnames(VECTOR_ELT(mmats, i)), 1);
292 :     SET_VECTOR_ELT(nms3, 0, cnms);
293 :     SET_VECTOR_ELT(nms3, 1, cnms);
294 :     SET_VECTOR_ELT(nms3, 2,
295 :     getAttrib(VECTOR_ELT(facs, i), R_LevelsSymbol));
296 :     dimnamesgets(VECTOR_ELT(bVar, i), duplicate(nms3));
297 :     }
298 :     SET_VECTOR_ELT(nms2, 0, xcols);
299 :     SET_VECTOR_ELT(nms2, 1, xcols);
300 :     dimnamesgets(GET_SLOT(x, Matrix_XtXSym), nms2);
301 :     dimnamesgets(GET_SLOT(x, Matrix_RXXSym), nms2);
302 :     UNPROTECT(2);
303 :     return R_NilValue;
304 :     }
305 :    
306 : bates 176 /**
307 :     * Update the numerical entries x, ZtX, and XtX in an ssclme object
308 :     * according to a set of model matrices.
309 :     *
310 :     * @param x pointer to an ssclme object
311 :     * @param facs pointer to a list of grouping factors
312 :     * @param mmats pointer to a list of model matrices
313 :     *
314 :     * @return NULL
315 :     */
316 : bates 22 SEXP
317 : bates 10 ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)
318 :     {
319 : bates 22 SEXP bVar = GET_SLOT(x, Matrix_bVarSym);
320 : bates 10 int
321 :     *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
322 :     *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
323 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
324 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
325 : bates 143 *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
326 : bates 10 i, j, k,
327 :     ione = 1,
328 :     nf = length(mmats) - 1,
329 :     nobs = nc[nf + 1],
330 :     nzcol = Gp[nf],
331 :     nz = Ap[nzcol],
332 :     pp1 = nc[nf];
333 :     double
334 :     **Z = Calloc(nf + 1, double *),
335 :     *Ax = REAL(GET_SLOT(x, Matrix_xSym)),
336 :     *XtX = REAL(GET_SLOT(x, Matrix_XtXSym)),
337 :     *ZtX = REAL(GET_SLOT(x, Matrix_ZtXSym)),
338 :     one = 1.0,
339 :     zero = 0.0;
340 :    
341 :     for (i = 0; i <= nf; i++) {
342 :     int *dims = INTEGER(getAttrib(VECTOR_ELT(mmats, i), R_DimSymbol)),
343 :     nci = nc[i];
344 :     if (nobs != dims[0])
345 :     error("Expected %d rows in the %d'th model matrix. Got %d",
346 :     nobs, i+1, dims[0]);
347 :     if (nci != dims[1])
348 :     error("Expected %d columns in the %d'th model matrix. Got %d",
349 :     nci, i+1, dims[1]);
350 :     Z[i] = REAL(VECTOR_ELT(mmats, i));
351 :     }
352 :     /* Create XtX - X is Z[nf] */
353 :     F77_CALL(dsyrk)("U", "T", nc+nf, &nobs, &one,
354 :     Z[nf], &nobs, &zero, XtX, nc + nf);
355 :     /* Zero the accumulators */
356 :     memset((void *) ZtX, 0, sizeof(double) * nzcol * pp1);
357 :     memset((void *) Ax, 0, sizeof(double) * nz);
358 :     for (j = 0; j < nf; j++) { /* Create ZtX */
359 :     int *fpj = INTEGER(VECTOR_ELT(facs, j)), ncj = nc[j],
360 :     Ncj = ncj > 1;
361 :     double
362 :     *bVj = REAL(VECTOR_ELT(bVar, j)),
363 :     *Zj = Z[j],
364 :     *zxj = ZtX + Gp[j];
365 :    
366 :     if (Ncj) { /* bVj will accumulate Z'Z blocks */
367 :     memset(bVj, 0, sizeof(double) * ncj * (Gp[j+1]-Gp[j]));
368 :     }
369 :     for (i = 0; i < nobs; i++) { /* accumulate diagonal of ZtZ */
370 :     int fpji = fpj[i] - 1, /* factor indices are 1-based */
371 :     dind = Ap[Gp[j] + fpji * ncj + 1] - 1;
372 :     if (Ai[dind] != (Gp[j] + fpji * ncj))
373 :     error("logic error in ssclme_update_mm");
374 :     if (Ncj) { /* use bVar to accumulate */
375 :     F77_CALL(dsyrk)("U", "T", &ncj, &ione, &one, Zj+i,
376 :     &nobs, &one, bVj + fpji*ncj*ncj, &ncj);
377 :     } else { /* update scalars directly */
378 :     Ax[dind] += Zj[i] * Zj[i];
379 :     }
380 :     /* update rows of Z'X */
381 :     F77_CALL(dgemm)("T", "N", &ncj, &pp1, &ione, &one,
382 :     Zj + i, &nobs, Z[nf] + i, &nobs,
383 :     &one, zxj + fpji * ncj, &nzcol);
384 :     }
385 :     if (Ncj) bVj_to_A(ncj, Gp[j], Gp[j+1], bVj, Ap, Ai, Ax);
386 :     for (k = j+1; k < nf; k++) { /* off-diagonals */
387 :     int *fpk = INTEGER(VECTOR_ELT(facs, k)),
388 :     *Apk = Ap + Gp[k],
389 : bates 111 nck = nc[k],
390 :     scalar = ncj == 1 && nck == 1;
391 : bates 10 double
392 : bates 245 *Zk = Z[k], *work = (double *) NULL;
393 :    
394 : bates 111 if (!scalar) work = Calloc(ncj * nck, double);
395 : bates 10 for (i = 0; i < nobs; i++) {
396 :     int ii, ind = -1, fpji = fpj[i] - 1,
397 :     row = Gp[j] + fpji * ncj,
398 :     fpki = fpk[i] - 1,
399 : bates 111 lastind = Apk[fpki*nck + 1];
400 :     for (ii = Apk[fpki*nck]; ii < lastind; ii++) {
401 : bates 10 if (Ai[ii] == row) {
402 :     ind = ii;
403 :     break;
404 :     }
405 :     }
406 :     if (ind < 0) error("logic error in ssclme_update_mm");
407 : bates 111 if (scalar) { /* update scalars directly */
408 : bates 108 Ax[ind] += Zj[i] * Zk[i];
409 : bates 111 } else {
410 :     int jj, offset = ind - Apk[fpki * nck];
411 :     F77_CALL(dgemm)("T", "N", &ncj, &nck, &ione, &one,
412 :     Zj + i, &nobs, Zk + i, &nobs,
413 :     &zero, work, &ncj);
414 :     for (jj = 0; jj < nck; jj++) {
415 :     ind = Apk[fpki * nck + jj] + offset;
416 :     if (Ai[ind] != row)
417 :     error("logic error in ssclme_update_mm");
418 :     for (ii = 0; ii < ncj; ii++) {
419 :     Ax[ind++] += work[jj * ncj + ii];
420 :     }
421 :     }
422 : bates 10 }
423 :     }
424 : bates 111 if (!scalar) Free(work);
425 : bates 10 }
426 :     }
427 :     Free(Z);
428 : bates 22 ssclme_transfer_dimnames(x, facs, mmats);
429 : bates 143 status[0] = status[1] = 0;
430 : bates 10 return R_NilValue;
431 :     }
432 :    
433 : bates 176 /**
434 :     * Inflate Z'Z according to Omega and create the factorization LDL'
435 :     *
436 :     * @param x pointer to an ssclme object
437 :     *
438 :     * @return NULL
439 :     */
440 : bates 97 SEXP ssclme_inflate_and_factor(SEXP x)
441 : bates 10 {
442 :     SEXP
443 : bates 97 GpSlot = GET_SLOT(x, Matrix_GpSym),
444 :     Omega = GET_SLOT(x, Matrix_OmegaSym);
445 :     int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
446 : bates 10 int
447 : bates 97 *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
448 :     *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
449 : bates 10 *Flag = Calloc(n, int),
450 :     *Gp = INTEGER(GpSlot),
451 :     *Lnz = Calloc(n, int),
452 :     *Pattern = Calloc(n, int),
453 : bates 97 *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
454 : bates 10 j,
455 :     nf = length(GpSlot) - 1;
456 :     double
457 : bates 97 *D = REAL(GET_SLOT(x, Matrix_DSym)),
458 :     *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
459 : bates 10 *Y = Calloc(n, double),
460 :     *xcp = Calloc(Ap[n], double);
461 :    
462 : bates 97 Memcpy(xcp, REAL(GET_SLOT(x, Matrix_xSym)), Ap[n]);
463 : bates 10 for (j = 0; j < nf; j++) {
464 :     int diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];
465 :     double *omgj = REAL(VECTOR_ELT(Omega, j));
466 :    
467 :     for (i = Gp[j]; i < G2; i += ncj) {
468 :     for (k = 0; k < ncj; k++) {
469 :     diag = Ap[i + k + 1] - 1;
470 :     if (Ai[diag] != i+k)
471 :     error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",
472 :     diag, i+k, Ai[diag]);
473 :     for (ii = 0; ii <= k; ii++) {
474 :     xcp[diag + ii - k] += omgj[k*ncj + ii];
475 :     }
476 :     }
477 :     }
478 :     }
479 : bates 372 j = R_ldl_numeric(n, Ap, Ai, xcp,
480 : bates 97 INTEGER(GET_SLOT(x, Matrix_LpSym)),
481 :     INTEGER(GET_SLOT(x, Matrix_ParentSym)),
482 : bates 372 INTEGER(GET_SLOT(x, Matrix_LiSym)),
483 : bates 97 REAL(GET_SLOT(x, Matrix_LxSym)),
484 : bates 372 D, (int *) NULL, (int *) NULL); /* P & Pinv */
485 : bates 10 if (j != n)
486 :     error("rank deficiency of ZtZ+W detected at column %d",
487 :     j + 1);
488 :     for (j = 0; j < n; j++) DIsqrt[j] = 1./sqrt(D[j]);
489 :     Free(Lnz); Free(Flag); Free(Pattern); Free(Y); Free(xcp);
490 :     return R_NilValue;
491 :     }
492 :    
493 : bates 176
494 :     /**
495 :     * If status[["factored"]] is FALSE, create and factor Z'Z+Omega, then
496 :     * create RZX and RXX, the deviance components, and the value of the
497 :     * deviance for both ML and REML.
498 :     *
499 :     * @param x pointer to an ssclme object
500 :     *
501 :     * @return NULL
502 :     */
503 : bates 97 SEXP ssclme_factor(SEXP x)
504 : bates 10 {
505 : bates 97 int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
506 : bates 10
507 :     if (!status[0]) {
508 :     SEXP
509 : bates 97 GpSlot = GET_SLOT(x, Matrix_GpSym),
510 :     Omega = GET_SLOT(x, Matrix_OmegaSym);
511 : bates 10 int
512 :     *Gp = INTEGER(GpSlot),
513 : bates 97 *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
514 :     *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
515 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
516 : bates 10 i,
517 : bates 97 n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1],
518 : bates 10 nf = length(GpSlot) - 1,
519 :     nobs = nc[nf + 1],
520 :     nreml = nobs + 1 - nc[nf],
521 :     pp1 = nc[nf],
522 :     pp2 = pp1 + 1;
523 :     double
524 : bates 97 *D = REAL(GET_SLOT(x, Matrix_DSym)),
525 :     *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
526 :     *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
527 :     *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
528 :     *RZX = REAL(GET_SLOT(x, Matrix_RZXSym)),
529 :     *dcmp = REAL(getAttrib(x, Matrix_devCompSym)),
530 :     *deviance = REAL(getAttrib(x, Matrix_devianceSym)),
531 : bates 10 minus1 = -1.,
532 :     one = 1.;
533 :    
534 : bates 97 ssclme_inflate_and_factor(x);
535 : bates 10 /* Accumulate logdet of ZtZ+W */
536 :     dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;
537 :     for (i = 0; i < n; i++) dcmp[0] += log(D[i]);
538 :     /* Accumulate logdet of W */
539 :     for (i = 0; i < nf; i++) {
540 :     int nci = nc[i],
541 :     mi = (Gp[i+1] - Gp[i])/nci;
542 :    
543 :     if (nci < 2) {
544 :     dcmp[1] += mi * log(REAL(VECTOR_ELT(Omega, i))[0]);
545 :     } else {
546 :     int j;
547 :     double
548 :     *tmp = Calloc(nci * nci, double),
549 :     accum = 0.;
550 :     F77_CALL(dpotrf)("U", &nci,
551 :     Memcpy(tmp, REAL(VECTOR_ELT(Omega, i)),
552 :     nci * nci),
553 :     &nci, &j);
554 :     if (j)
555 :     error("Omega[%d] is not positive definite", i + 1);
556 :     for (j = 0; j < nci; j++) {
557 :     accum += 2 * log(tmp[j * (nci + 1)]);
558 :     }
559 :     dcmp[1] += mi * accum;
560 :     Free(tmp);
561 :     }
562 :     }
563 :     /* ldl_lsolve on Z'X */
564 : bates 97 Memcpy(RZX, REAL(GET_SLOT(x, Matrix_ZtXSym)), n * pp1);
565 : bates 10 for (i = 0; i < pp1; i++) {
566 :     int j;
567 :     double *RZXi = RZX + i * n;
568 : bates 372 R_ldl_lsolve(n, RZXi, Lp, Li, Lx);
569 : bates 10 for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
570 :     }
571 :     /* downdate and factor X'X */
572 : bates 97 Memcpy(RXX, REAL(GET_SLOT(x, Matrix_XtXSym)), pp1 * pp1);
573 : bates 10 F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
574 :     RZX, &n, &one, RXX, &pp1);
575 :     F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
576 : bates 122 if (i) {
577 :     warning("Could not factor downdated X'X, code %d", i);
578 :     dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;
579 :     } else {
580 : bates 10 /* logdet of RXX */
581 : bates 122 for (i = 0; i < (pp1 - 1); i++)
582 :     dcmp[2] += 2 * log(RXX[i*pp2]);
583 : bates 10 /* logdet of Ryy */
584 : bates 122 dcmp[3] = 2. * log(RXX[pp1 * pp1 - 1]);
585 :     deviance[0] = /* ML criterion */
586 :     dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));
587 :     deviance[1] = dcmp[0] - dcmp[1] + /* REML */
588 :     dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));
589 :     }
590 : bates 10 status[0] = 1;
591 :     status[1] = 0;
592 :     }
593 :     return R_NilValue;
594 :     }
595 :    
596 : bates 176 /**
597 :     * Return the position of probe in the sorted index vector ind. It is
598 :     * known that the position is greater than or equal to start so a linear
599 :     * search from start is used.
600 :     *
601 :     * @param probe value to be matched
602 :     * @param start index at which to start
603 :     * @param ind vector of indices
604 :     *
605 :     * @return index of the entry matching probe
606 :     */
607 : bates 10 static
608 :     int ldl_update_ind(int probe, int start, const int ind[])
609 :     {
610 :     while (ind[start] < probe) start++;
611 :     if (ind[start] > probe) error("logic error in ldl_inverse");
612 :     return start;
613 :     }
614 :    
615 :     /**
616 : bates 176 * Update the diagonal blocks of the inverse of LDL' (=Z'Z+W). The
617 :     * lower Cholesky factors of the updated blocks are stored in the bVar
618 :     * slot.
619 : bates 10 *
620 :     * @param x pointer to an ssclme object
621 :     *
622 :     * @return R_NilValue (x is updated in place)
623 :    
624 :     */
625 : bates 97 static
626 : bates 245 void ldl_inverse(SEXP x)
627 : bates 10 {
628 :     SEXP
629 :     Gpsl = GET_SLOT(x, Matrix_GpSym),
630 :     bVar = GET_SLOT(x, Matrix_bVarSym);
631 :     int *Gp = INTEGER(Gpsl),
632 : bates 159 *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym)),
633 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
634 : bates 10 i,
635 :     nf = length(Gpsl) - 1,
636 : bates 159 nzc = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
637 :     int maxod = Parent[nzc];
638 :     double *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym));
639 : bates 10
640 :     ssclme_factor(x);
641 : bates 159 if (maxod == 0) { /* L and L^{-1} are the identity */
642 : bates 10 for (i = 0; i < nf; i++) {
643 :     Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],
644 :     Gp[i+1] - Gp[i]);
645 :     }
646 : bates 159 } else {
647 :     int *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
648 :     *Li = INTEGER(GET_SLOT(x, Matrix_LiSym));
649 :    
650 :     double one = 1.0, zero = 0.,
651 :     *Lx = REAL(GET_SLOT(x, Matrix_LxSym));
652 : bates 157
653 : bates 10 for (i = 0; i < nf; i++) {
654 : bates 158 int j, jj, k, kk, nci = nc[i], nr, p, p2, pj, pp,
655 : bates 157 m = maxod + 1,
656 : bates 159 *ind = Calloc(m, int), G1 = Gp[i], G2 = Gp[i+1];
657 : bates 10 double
658 :     *tmp = Calloc(m * nci, double),
659 : bates 159 *bVi = REAL(VECTOR_ELT(bVar, i));
660 : bates 10
661 : bates 159 /* initialize bVi to zero */
662 :     memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
663 :    
664 :     for (j = G1; j < G2; j += nci) {
665 : bates 158 kk = 0; /* ind gets indices of non-zeros */
666 : bates 10 jj = j; /* in this block of columns */
667 :     while (jj >= 0) {
668 :     ind[kk++] = jj;
669 :     jj = Parent[jj];
670 :     }
671 :     nr = kk; /* number of non-zeros in this block */
672 :     while (kk < m) ind[kk++] = nzc; /* placeholders */
673 :    
674 :     for (k = 0; k < nci; k++) {
675 :     double *ccol = tmp + k * nr;
676 :    
677 : bates 157 for (kk = 0; kk < nr; kk++) ccol[kk] = 0.;
678 :     ccol[k] = 1.; /* initialize from unit diagonal */
679 : bates 10 for (jj = j + k; jj >= 0; jj = Parent[jj]) {
680 :     p2 = Lp[jj+1];
681 : bates 158 pp = pj = ldl_update_ind(jj, 0, ind);
682 : bates 10 for (p = Lp[jj]; p < p2; p++) {
683 :     pp = ldl_update_ind(Li[p], pp, ind);
684 : bates 158 ccol[pp] -= Lx[p] * ccol[pj];
685 : bates 10 }
686 :     }
687 :     }
688 : bates 157
689 :     for (kk = 0; kk < nr; kk++) { /* scale rows */
690 : bates 10 for (k = 0; k < nci; k++) {
691 :     tmp[k * nr + kk] *= DIsqrt[ind[kk]];
692 :     }
693 :     }
694 : bates 113 F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,
695 : bates 159 &zero, bVi + (j - G1)*nci, &nci);
696 :     F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,
697 :     &nci, &jj);
698 :     if (jj) /* should never happen */
699 : bates 10 error(
700 : bates 159 "Rank deficient variance matrix at group %d, level %d, error code %d",
701 :     i + 1, j + 1, jj);
702 : bates 10 }
703 :     Free(tmp); Free(ind);
704 :     }
705 :     }
706 :     }
707 : bates 101
708 : bates 176 /**
709 :     * If necessary, factor Z'Z+Omega, ZtX, and XtX then, if necessary,
710 :     * form RZX, RXX, and bVar for the inverse of the Cholesky factor.
711 :     *
712 :     * @param x pointer to an ssclme object
713 :     *
714 :     * @return NULL (x is updated in place)
715 :     */
716 : bates 97 SEXP ssclme_invert(SEXP x)
717 : bates 10 {
718 : bates 97 int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
719 :     if (!status[0]) ssclme_factor(x);
720 : bates 122 if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))
721 :     error("Unable to invert singular factor of downdated X'X");
722 : bates 10 if (!status[1]) {
723 :     SEXP
724 : bates 97 RZXsl = GET_SLOT(x, Matrix_RZXSym);
725 : bates 10 int
726 :     *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
727 : bates 97 *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
728 :     *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
729 : bates 10 i,
730 :     n = dims[0],
731 :     pp1 = dims[1];
732 :     double
733 : bates 97 *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
734 :     *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
735 :     *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
736 : bates 10 *RZX = REAL(RZXsl),
737 :     one = 1.;
738 :    
739 :     F77_CALL(dtrtri)("U", "N", &pp1, RXX, &pp1, &i);
740 :     if (i)
741 :     error("DTRTRI returned error code %d", i);
742 :     F77_CALL(dtrmm)("R", "U", "N", "N", &n, &pp1, &one,
743 :     RXX, &pp1, RZX, &n);
744 :     for (i = 0; i < pp1; i++) {
745 :     int j; double *RZXi = RZX + i * n;
746 :     for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
747 : bates 372 R_ldl_ltsolve(n, RZXi, Lp, Li, Lx);
748 : bates 10 }
749 : bates 97 ldl_inverse(x);
750 : bates 10 status[1] = 1;
751 :     }
752 :     return R_NilValue;
753 :     }
754 :    
755 : bates 176 /**
756 :     * Create and insert initial values for Omega_i.
757 :     *
758 :     * @param x pointer to an ssclme object
759 :     *
760 :     * @return NULL
761 :     */
762 : bates 10 SEXP ssclme_initial(SEXP x)
763 :     {
764 :     SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),
765 :     Omg = GET_SLOT(x, Matrix_OmegaSym);
766 :     int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
767 :     *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
768 :     *Gp = INTEGER(Gpsl),
769 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
770 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
771 :     i, nf = length(Gpsl) - 1;
772 :    
773 :     double *Ax = REAL(GET_SLOT(x, Matrix_xSym));
774 :    
775 :     for (i = 0; i < nf; i++) {
776 :     int
777 :     Gpi = Gp[i],
778 :     j, k,
779 :     nci = nc[i],
780 :     ncip1 = nci + 1,
781 :     p2 = Gp[i+1];
782 :     double
783 :     mi = 0.375 * ((double) nci)/((double) (p2 - Gpi)),
784 :     *mm = REAL(VECTOR_ELT(Omg, i));
785 :    
786 :     memset((void *) mm, 0, sizeof(double) * nci * nci);
787 :     for (j = Gpi; j < p2; j += nci) {
788 :     for (k = 0; k < nci; k++) {
789 :     int jk = j+k, jj = Ap[jk+1] - 1;
790 :     if (Ai[jj] != jk) error("malformed ZtZ structure");
791 :     mm[k * ncip1] += Ax[jj] * mi;
792 :     }
793 :     }
794 :     }
795 :     status[0] = status[1] = 0;
796 :     return R_NilValue;
797 :     }
798 :    
799 :     /**
800 :     * Extract the conditional estimates of the fixed effects
801 :     *
802 :     * @param x Pointer to an ssclme object
803 :     *
804 :     * @return a numeric vector containing the conditional estimates of
805 :     * the fixed effects
806 :     */
807 :     SEXP ssclme_fixef(SEXP x)
808 :     {
809 :     SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
810 :     int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1];
811 :     int j, p = pp1 - 1;
812 :     SEXP val = PROTECT(allocVector(REALSXP, p));
813 :     double
814 :     *RXX = REAL(RXXsl),
815 :     *beta = REAL(val),
816 :     nryyinv; /* negative ryy-inverse */
817 :    
818 :     ssclme_invert(x);
819 :     Memcpy(beta, RXX + p * pp1, p);
820 :     nryyinv = -RXX[pp1*pp1 - 1];
821 :     for (j = 0; j < p; j++) beta[j] /= nryyinv;
822 :     UNPROTECT(1);
823 :     return val;
824 :     }
825 :    
826 :     /**
827 :     * Extract the conditional modes of the random effects.
828 :     *
829 :     * @param x Pointer to an ssclme object
830 :     *
831 :     * @return a vector containing the conditional modes of the random effects
832 :     */
833 :     SEXP ssclme_ranef(SEXP x)
834 :     {
835 : bates 21 SEXP RZXsl = GET_SLOT(x, Matrix_RZXSym),
836 :     GpSl = GET_SLOT(x, Matrix_GpSym);
837 : bates 10 int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
838 : bates 21 *Gp = INTEGER(GpSl),
839 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
840 :     i, j,
841 : bates 10 n = dims[0],
842 : bates 21 nf = length(GpSl) - 1,
843 : bates 10 pp1 = dims[1],
844 :     p = pp1 - 1;
845 : bates 21 SEXP val = PROTECT(allocVector(VECSXP, nf));
846 : bates 10 double
847 : bates 21 *b = REAL(RZXsl) + n * p,
848 : bates 10 ryyinv; /* ryy-inverse */
849 :    
850 :     ssclme_invert(x);
851 :     ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];
852 : bates 21 for (i = 0; i < nf; i++) {
853 : bates 34 int nci = nc[i], Mi = Gp[i+1] - Gp[i];
854 : bates 21 double *mm;
855 :    
856 : bates 34 SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, Mi/nci));
857 : bates 21 mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);
858 : bates 34 b += Mi;
859 : bates 21 for (j = 0; j < Mi; j++) mm[j] /= ryyinv;
860 :     }
861 : bates 10 UNPROTECT(1);
862 :     return val;
863 :     }
864 : bates 34
865 : bates 10 /**
866 :     * Extract the ML or REML conditional estimate of sigma
867 :     *
868 :     * @param x pointer to an ssclme object
869 :     * @param REML logical scalar - TRUE if REML estimates are requested
870 :     *
871 :     * @return numeric scalar
872 :     */
873 :     SEXP ssclme_sigma(SEXP x, SEXP REML)
874 :     {
875 :     SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
876 :     int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],
877 : bates 122 nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))[
878 :     length(GET_SLOT(x, Matrix_GpSym))];
879 : bates 10
880 :     ssclme_invert(x);
881 :     return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *
882 :     sqrt((double)(asLogical(REML) ?
883 :     nobs + 1 - pp1 : nobs))));
884 :     }
885 :    
886 : bates 176 /**
887 :     * Calculate the length of the parameter vector, which is called coef
888 :     * for historical reasons.
889 :     *
890 :     * @param nf number of factors
891 :     * @param nc number of columns in the model matrices for each factor
892 :     *
893 :     * @return total length of the coefficient vector
894 :     */
895 : bates 10 static
896 :     int coef_length(int nf, const int nc[])
897 :     {
898 :     int i, ans = 0;
899 :     for (i = 0; i < nf; i++) ans += (nc[i] * (nc[i] + 1))/2;
900 :     return ans;
901 :     }
902 :    
903 :     /**
904 : bates 176 * Extract the upper triangles of the Omega matrices. These aren't
905 :     * "coefficients" but the extractor is called coef for historical
906 :     * reasons. Within each group these values are in the order of the
907 :     * diagonal entries first then the strict upper triangle in row
908 :     * order.
909 : bates 10 *
910 :     * @param x pointer to an ssclme object
911 :     *
912 :     * @return numeric vector of the values in the upper triangles of the
913 :     * Omega matrices
914 :     */
915 : bates 247 SEXP ssclme_coef(SEXP x, SEXP Unconstr)
916 : bates 10 {
917 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
918 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
919 : bates 247 i, nf = length(Omega), unc = asLogical(Unconstr), vind;
920 : bates 10 SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
921 :     double *vv = REAL(val);
922 :    
923 : bates 247 vind = 0; /* index in vv */
924 : bates 10 for (i = 0; i < nf; i++) {
925 : bates 247 int nci = nc[i], ncip1 = nci + 1;
926 : bates 104 if (nci == 1) {
927 : bates 247 vv[vind++] = (unc ?
928 :     log(REAL(VECTOR_ELT(Omega, i))[0]) :
929 :     REAL(VECTOR_ELT(Omega, i))[0]);
930 : bates 104 } else {
931 : bates 247 if (unc) { /* L log(D) L' factor of Omega[i,,] */
932 :     int j, k, ncisq = nci * nci;
933 :     double *tmp = Memcpy(Calloc(ncisq, double),
934 :     REAL(VECTOR_ELT(Omega, i)), ncisq);
935 :     F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
936 :     if (j) /* should never happen */
937 :     error("DPOTRF returned error code %d on Omega[[%d]]",
938 :     j, i+1);
939 :     for (j = 0; j < nci; j++) {
940 :     double diagj = tmp[j * ncip1];
941 :     vv[vind++] = 2. * log(diagj);
942 :     for (k = j + 1; k < nci; k++) {
943 :     tmp[j + k * nci] /= diagj;
944 :     }
945 : bates 104 }
946 : bates 247 for (j = 0; j < nci; j++) {
947 :     for (k = j + 1; k < nci; k++) {
948 :     vv[vind++] = tmp[j + k * nci];
949 :     }
950 :     }
951 :     Free(tmp);
952 :     } else { /* upper triangle of Omega[i,,] */
953 :     int j, k, odind = vind + nci;
954 :     double *omgi = REAL(VECTOR_ELT(Omega, i));
955 :    
956 :     for (j = 0; j < nci; j++) {
957 :     vv[vind++] = omgi[j * ncip1];
958 :     for (k = j + 1; k < nci; k++) {
959 :     vv[odind++] = omgi[k*nci + j];
960 :     }
961 :     }
962 :     vind = odind;
963 : bates 10 }
964 :     }
965 :     }
966 :     UNPROTECT(1);
967 :     return val;
968 :     }
969 :    
970 :     /**
971 : bates 104 * Extract the unconstrained parameters that determine the
972 : bates 176 * Omega matrices. (Called coef for historical reasons.) The
973 :     * unconstrained parameters are derived from the LDL' decomposition of
974 :     * Omega_i. The first nc[i] entries in each group are the diagonals
975 :     * of log(D) followed by the strict lower triangle of L in column
976 :     * order.
977 : bates 19 *
978 :     * @param x pointer to an ssclme object
979 :     *
980 : bates 104 * @return numeric vector of unconstrained parameters that determine the
981 : bates 19 * Omega matrices
982 :     */
983 :     SEXP ssclme_coefUnc(SEXP x)
984 :     {
985 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
986 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
987 :     i, nf = length(Omega), vind;
988 :     SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
989 :     double *vv = REAL(val);
990 :    
991 :     vind = 0;
992 :     for (i = 0; i < nf; i++) {
993 :     int nci = nc[i];
994 :     if (nci == 1) {
995 :     vv[vind++] = log(REAL(VECTOR_ELT(Omega, i))[0]);
996 :     } else {
997 :     int j, k, ncip1 = nci + 1, ncisq = nci * nci;
998 :     double *tmp = Memcpy(Calloc(ncisq, double),
999 :     REAL(VECTOR_ELT(Omega, i)), ncisq);
1000 :     F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1001 :     if (j) /* should never happen */
1002 : bates 104 error("DPOTRF returned error code %d on Omega[[%d]]",
1003 :     j, i+1);
1004 : bates 19 for (j = 0; j < nci; j++) {
1005 :     double diagj = tmp[j * ncip1];
1006 :     vv[vind++] = 2. * log(diagj);
1007 :     for (k = j + 1; k < nci; k++) {
1008 :     tmp[j + k * nci] /= diagj;
1009 :     }
1010 :     }
1011 :     for (j = 0; j < nci; j++) {
1012 :     for (k = j + 1; k < nci; k++) {
1013 :     vv[vind++] = tmp[j + k * nci];
1014 :     }
1015 :     }
1016 :     Free(tmp);
1017 :     }
1018 :     }
1019 :     UNPROTECT(1);
1020 :     return val;
1021 :     }
1022 :    
1023 :     /**
1024 : bates 104 * Assign the Omega matrices from the unconstrained parameterization.
1025 : bates 19 *
1026 :     * @param x pointer to an ssclme object
1027 :     * @param coef pointer to an numeric vector of appropriate length
1028 :     *
1029 :     * @return R_NilValue
1030 :     */
1031 :     SEXP ssclme_coefGetsUnc(SEXP x, SEXP coef)
1032 :     {
1033 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1034 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1035 :     cind, i, nf = length(Omega),
1036 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1037 :     double *cc = REAL(coef);
1038 :    
1039 :     if (length(coef) != coef_length(nf, nc) || !isReal(coef))
1040 :     error("coef must be a numeric vector of length %d",
1041 :     coef_length(nf, nc));
1042 :     cind = 0;
1043 :     for (i = 0; i < nf; i++) {
1044 :     int nci = nc[i];
1045 :     if (nci == 1) {
1046 :     REAL(VECTOR_ELT(Omega, i))[0] = exp(cc[cind++]);
1047 :     } else {
1048 :     int odind = cind + nci, /* off-diagonal index */
1049 :     j, k,
1050 :     ncip1 = nci + 1,
1051 :     ncisq = nci * nci;
1052 :     double
1053 :     *omgi = REAL(VECTOR_ELT(Omega, i)),
1054 :     *tmp = Calloc(ncisq, double),
1055 : bates 104 diagj, one = 1., zero = 0.;
1056 :    
1057 : bates 19 memset(omgi, 0, sizeof(double) * ncisq);
1058 :     for (j = 0; j < nci; j++) {
1059 : bates 104 tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1060 : bates 19 for (k = j + 1; k < nci; k++) {
1061 : bates 104 tmp[k*nci + j] = cc[odind++] * diagj;
1062 : bates 19 }
1063 :     }
1064 : bates 104 F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1065 :     tmp, &nci, &zero, omgi, &nci);
1066 : bates 19 Free(tmp);
1067 :     cind = odind;
1068 :     }
1069 :     }
1070 :     status[0] = status[1] = 0;
1071 :     return x;
1072 :     }
1073 :    
1074 :     /**
1075 : bates 10 * Assign the upper triangles of the Omega matrices.
1076 : bates 104 * (Called coef for historical reasons.)
1077 : bates 10 *
1078 :     * @param x pointer to an ssclme object
1079 :     * @param coef pointer to an numeric vector of appropriate length
1080 :     *
1081 :     * @return R_NilValue
1082 :     */
1083 : bates 247 SEXP ssclme_coefGets(SEXP x, SEXP coef, SEXP Unc)
1084 : bates 10 {
1085 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1086 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1087 : bates 247 *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1088 : bates 10 cind, i, nf = length(Omega),
1089 : bates 247 unc = asLogical(Unc);
1090 : bates 10 double *cc = REAL(coef);
1091 :    
1092 :     if (length(coef) != coef_length(nf, nc) || !isReal(coef))
1093 :     error("coef must be a numeric vector of length %d",
1094 :     coef_length(nf, nc));
1095 :     cind = 0;
1096 :     for (i = 0; i < nf; i++) {
1097 : bates 104 int nci = nc[i];
1098 :     if (nci == 1) {
1099 : bates 247 REAL(VECTOR_ELT(Omega, i))[0] = (unc ?
1100 :     exp(cc[cind++]) :
1101 :     cc[cind++]);
1102 : bates 104 } else {
1103 : bates 247 int odind = cind + nci, /* off-diagonal index */
1104 :     j, k,
1105 :     ncip1 = nci + 1,
1106 :     ncisq = nci * nci;
1107 :     double
1108 :     *omgi = REAL(VECTOR_ELT(Omega, i));
1109 :     if (unc) {
1110 :     double
1111 :     *tmp = Calloc(ncisq, double),
1112 :     diagj, one = 1., zero = 0.;
1113 :    
1114 :     memset(omgi, 0, sizeof(double) * ncisq);
1115 :     for (j = 0; j < nci; j++) {
1116 :     tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1117 :     for (k = j + 1; k < nci; k++) {
1118 :     tmp[k*nci + j] = cc[odind++] * diagj;
1119 :     }
1120 : bates 104 }
1121 : bates 247 F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1122 :     tmp, &nci, &zero, omgi, &nci);
1123 :     Free(tmp);
1124 :     } else {
1125 :     for (j = 0; j < nci; j++) {
1126 :     omgi[j * ncip1] = cc[cind++];
1127 :     for (k = j + 1; k < nci; k++) {
1128 :     omgi[k*nci + j] = cc[odind++];
1129 :     }
1130 :     }
1131 : bates 10 }
1132 : bates 104 cind = odind;
1133 : bates 10 }
1134 :     }
1135 :     status[0] = status[1] = 0;
1136 : bates 11 return x;
1137 : bates 10 }
1138 :    
1139 : bates 245
1140 : bates 176 /**
1141 : bates 245 * Returns the inverse of the updated Omega matrices for an ECME
1142 :     * iteration. These matrices are also used in the gradient calculation.
1143 :     *
1144 :     * @param x pointer to an ssclme object
1145 :     * @param REML indicator of REML being used
1146 :     * @param val pointer to a list of symmetric q_i by q_i matrices
1147 :     */
1148 :     static void
1149 :     common_ECME_gradient(SEXP x, int REML, SEXP val)
1150 :     {
1151 :     SEXP
1152 :     RZXsl = GET_SLOT(x, Matrix_RZXSym),
1153 :     bVar = GET_SLOT(x, Matrix_bVarSym);
1154 :     int
1155 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1156 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1157 :     i, j, n = *INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1158 :     nf = length(val), nobs = nc[nf + 1], p = nc[nf] - 1;
1159 :     double
1160 :     *RZX = REAL(RZXsl),
1161 :     *b = RZX + p * INTEGER(getAttrib(RZXsl, R_DimSymbol))[0],
1162 :     one = 1.0, zero = 0.0;
1163 :    
1164 :     ssclme_invert(x);
1165 :     for (i = 0; i < nf; i++) {
1166 :     int ki = Gp[i+1] - Gp[i],
1167 :     nci = nc[i],
1168 :     mi = ki/nci;
1169 :     double
1170 :     *vali = REAL(VECTOR_ELT(val, i)),
1171 :     alpha = (double)(REML ? (nobs-p) : nobs)/((double) mi);
1172 :    
1173 :     F77_CALL(dsyrk)("U", "N", &nci, &mi,
1174 :     &alpha, b + Gp[i], &nci,
1175 :     &zero, vali, &nci);
1176 :     alpha = 1./((double) mi);
1177 :     F77_CALL(dsyrk)("U", "N", &nci, &ki,
1178 :     &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1179 :     &one, vali, &nci);
1180 :     if (REML) {
1181 :     for (j = 0; j < p; j++) {
1182 :     F77_CALL(dsyrk)("U", "N", &nci, &mi,
1183 :     &alpha, RZX + Gp[i] + j*n, &nci,
1184 :     &one, vali, &nci);
1185 :     }
1186 :     }
1187 :     }
1188 :     }
1189 :    
1190 :     /**
1191 : bates 247 * Print the verbose output in the ECME iterations
1192 : bates 176 *
1193 :     * @param x pointer to an ssclme object
1194 : bates 247 * @param iter iteration number
1195 :     * @param REML indicator of whether REML is being used
1196 :     */
1197 :     static void EMsteps_verbose_print(SEXP x, int iter, int REML)
1198 :     {
1199 :     SEXP coef = PROTECT(ssclme_coef(x, ScalarLogical(0)));
1200 : bates 250 int i, lc = length(coef);
1201 :     double *cc = REAL(coef), *dev = REAL(GET_SLOT(x, Matrix_devianceSym));
1202 :    
1203 : bates 247
1204 :     ssclme_factor(x);
1205 :     if (iter == 0) Rprintf(" EM iterations\n");
1206 :     Rprintf("%3d %.3f", iter, dev[REML ? 1 : 0]);
1207 :     for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1208 :     Rprintf("\n");
1209 :     UNPROTECT(1);
1210 :     }
1211 :    
1212 :     /**
1213 :     * Perform ECME steps for the REML or ML criterion.
1214 :     *
1215 :     * @param x pointer to an ssclme object
1216 :     * @param nsteps pointer to an integer scalar - the number of ECME steps to perform
1217 : bates 176 * @param REMLp pointer to a logical scalar indicating if REML is to be used
1218 :     * @param verb pointer to a logical scalar indicating verbose mode
1219 :     *
1220 :     * @return NULL
1221 :     */
1222 : bates 11 SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1223 : bates 10 {
1224 :     SEXP
1225 : bates 245 Omega = GET_SLOT(x, Matrix_OmegaSym);
1226 : bates 10 int
1227 : bates 245 *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1228 : bates 10 *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1229 :     REML = asLogical(REMLp),
1230 :     i, info, iter,
1231 :     nEM = asInteger(nsteps),
1232 : bates 245 nf = length(Omega),
1233 : bates 10 verbose = asLogical(verb);
1234 :    
1235 : bates 247 if (verbose) EMsteps_verbose_print(x, 0, REML);
1236 : bates 166 for (iter = 0; iter < nEM; iter++) {
1237 : bates 245 common_ECME_gradient(x, REML, Omega);
1238 :     status[0] = status[1] = 0;
1239 : bates 10 for (i = 0; i < nf; i++) {
1240 : bates 245 double *vali = REAL(VECTOR_ELT(Omega, i));
1241 :    
1242 :     F77_CALL(dpotrf)("U", nc + i, vali, nc + i, &info);
1243 : bates 114 if (info)
1244 :     error("DPOTRF returned error code %d in Omega[[%d]] update",
1245 :     info, i + 1);
1246 : bates 245 F77_CALL(dpotri)("U", nc + i, vali, nc + i, &info);
1247 : bates 114 if (info)
1248 :     error("DPOTRI returned error code %d in Omega[[%d]] update",
1249 :     info, i + 1);
1250 : bates 10 }
1251 : bates 247 if (verbose) EMsteps_verbose_print(x, iter + 1, REML);
1252 : bates 10 }
1253 :     ssclme_factor(x);
1254 :     return R_NilValue;
1255 :     }
1256 : bates 11
1257 : bates 176 /**
1258 : bates 245 * Evaluate the gradient with respect to the indicators of the
1259 :     * positions in the Omega matrices.
1260 :     *
1261 :     * @param x Pointer to an ssclme object
1262 :     * @param REML indicator of whether REML is to be used
1263 :     * @param val Pointer to an object of the same structure as Omega
1264 :     */
1265 :     static void indicator_gradient(SEXP x, int REML, SEXP val)
1266 :     {
1267 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1268 :     int
1269 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1270 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1271 :     i, nf = length(Omega);
1272 :    
1273 :     common_ECME_gradient(x, REML, val);
1274 :     for (i = 0; i < nf; i++) {
1275 :     int info, nci = nc[i], mi = (Gp[i+1] - Gp[i])/nci;
1276 :     double
1277 :     *work = Memcpy((double *) Calloc(nci * nci, double),
1278 :     REAL(VECTOR_ELT(Omega, i)), nci * nci),
1279 :     alpha = (double) -mi, beta = (double) mi;
1280 :    
1281 :     F77_CALL(dpotrf)("U", &nci, work, &nci, &info);
1282 :     if (info)
1283 :     error("DPOTRF gave error code %d on Omega[[%d]]", info, i + 1);
1284 :     F77_CALL(dtrtri)("U", "N", &nci, work, &nci, &info);
1285 :     if (info)
1286 :     error("DPOTRI gave error code %d on Omega[[%d]]", info, i + 1);
1287 :     F77_CALL(dsyrk)("U", "N", &nci, &nci, &alpha, work, &nci,
1288 :     &beta, REAL(VECTOR_ELT(val, i)), &nci);
1289 :     Free(work);
1290 :     }
1291 :     }
1292 :    
1293 :     /**
1294 :     * Overwrite a gradient with respect to positions in Omega[[i]] by the
1295 :     * gradient with respect to the unconstrained parameters.
1296 :     *
1297 : bates 247 * @param grad pointer to a gradient wrt positions. Contents are overwritten.
1298 :     * @param Omega pointer to a list of symmetric matrices (upper storage).
1299 : bates 245 */
1300 :     static void unconstrained_gradient(SEXP grad, SEXP Omega)
1301 :     {
1302 :     int i, ii, j, nf = length(Omega);
1303 :     double one = 1.0, zero = 0.0;
1304 :    
1305 :     for (i = 0; i < nf; i++) {
1306 :     SEXP Omegai = VECTOR_ELT(Omega, i);
1307 :     int nci = *INTEGER(getAttrib(Omegai, R_DimSymbol)),
1308 :     ncisq = nci * nci, ncip1 = nci + 1;
1309 :     double *chol =
1310 :     Memcpy(Calloc(ncisq, double), REAL(Omegai), ncisq),
1311 :     *ai = REAL(VECTOR_ELT(grad, i)),
1312 :     *tmp = Calloc(ncisq, double);
1313 :    
1314 :     F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1315 :     if (j)
1316 :     error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1317 :     /* calculate and store grad[i,,] %*% t(R) */
1318 :     for (j = 0; j < nci; j++) {
1319 :     F77_CALL(dsymv)("U", &nci, &one, ai, &nci, chol + j, &nci,
1320 :     &zero, tmp + j, &nci);
1321 :     }
1322 :     /* full symmetric product gives diagonals */
1323 :     F77_CALL(dtrmm)("R", "U", "T", "N", &nci, &nci, &one, chol, &nci,
1324 :     Memcpy(ai, tmp, ncisq), &nci);
1325 : bates 247 /* overwrite upper triangle with gradients for positions in L' */
1326 : bates 245 for (ii = 1; ii < nci; ii++) {
1327 :     for (j = 0; j < ii; j++) {
1328 : bates 247 ai[j + ii*nci] = 2. * chol[j*ncip1] * tmp[j + ii*nci];
1329 :     ai[ii + j*nci] = 0.;
1330 : bates 245 }
1331 :     }
1332 :     Free(chol); Free(tmp);
1333 :     }
1334 :     }
1335 :    
1336 : bates 247 /**
1337 :     * Fills cvec with unlist(lapply(mList, function(el) el[upper.tri(el, strict = FALSE)]))
1338 :     *
1339 :     * @param mList pointer to a list of REAL matrices
1340 :     * @param nc number of columns in each matrix
1341 :     * @param cvec pointer to REAL vector to receive the result
1342 :     */
1343 :     static void upperTriList_to_vector(SEXP mList, int *nc, SEXP cvec)
1344 : bates 245 {
1345 : bates 247 int i, ii, j, nf = length(mList), pos = 0;
1346 : bates 245
1347 : bates 247 pos = 0; /* position in vector */
1348 :     for (i = 0; i < nf; i++) {
1349 :     SEXP omgi = VECTOR_ELT(mList, i);
1350 :     int nci = nc[i];
1351 :     for (j = 0; j < nci; j++) {
1352 :     for (ii = 0; ii <= j; ii++) {
1353 :     REAL(cvec)[pos++] = REAL(omgi)[ii + j * nci];
1354 :     }
1355 :     }
1356 :     }
1357 : bates 245 }
1358 : bates 250
1359 : bates 245 /**
1360 : bates 176 * Return the gradient of the ML or REML deviance.
1361 :     *
1362 :     * @param x pointer to an ssclme object
1363 :     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1364 :     * @param Uncp pointer to a logical scalar indicating if the unconstrained parameterization is to be used
1365 : bates 247 * @param Uncp pointer to a logical scalar indicating if result should be a single vector
1366 : bates 176 *
1367 : bates 247 * @return pointer to the gradient as a list of matrices or as a vector.
1368 : bates 176 */
1369 : bates 247 SEXP ssclme_grad(SEXP x, SEXP REMLp, SEXP Unc, SEXP OneVector)
1370 :     {
1371 :     SEXP ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1372 :    
1373 :     indicator_gradient(x, asLogical(REMLp), ans);
1374 :     if (asLogical(Unc))
1375 :     unconstrained_gradient(ans, GET_SLOT(x, Matrix_OmegaSym));
1376 :     if (asLogical(OneVector)) {
1377 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));
1378 :     SEXP val = PROTECT(allocVector(REALSXP, coef_length(length(ans), nc)));
1379 :    
1380 : bates 250 upperTriList_to_vector(ans, nc, val);
1381 : bates 247 UNPROTECT(2);
1382 :     return val;
1383 :     }
1384 :     UNPROTECT(1);
1385 :     return ans;
1386 :     }
1387 :    
1388 : bates 101 SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1389 : bates 100 {
1390 :     SEXP
1391 : bates 101 Omega = GET_SLOT(x, Matrix_OmegaSym),
1392 : bates 100 RZXsl = GET_SLOT(x, Matrix_RZXSym),
1393 :     bVar = GET_SLOT(x, Matrix_bVarSym);
1394 :     int
1395 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1396 : bates 176 *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1397 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1398 : bates 100 REML = asLogical(REMLp),
1399 : bates 104 cind, i, n = dims[0],
1400 : bates 105 nf = length(Omega),
1401 : bates 113 nobs, p, pp1 = dims[1],
1402 : bates 101 uncst = asLogical(Uncp);
1403 : bates 100 double
1404 :     *RZX = REAL(RZXsl),
1405 :     *b,
1406 :     alpha,
1407 :     one = 1.;
1408 : bates 105 SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1409 : bates 100
1410 : bates 125 ssclme_factor(x);
1411 :     if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {
1412 :     int ncoef = coef_length(nf, nc);
1413 :     for (i = 0; i < ncoef; i++) REAL(ans)[i] = NA_REAL;
1414 :     UNPROTECT(1);
1415 :     return ans;
1416 :     }
1417 : bates 104 nobs = nc[nf + 1];
1418 : bates 100 p = pp1 - 1;
1419 :     b = RZX + p * n;
1420 :     ssclme_invert(x);
1421 : bates 104 cind = 0;
1422 : bates 100 for (i = 0; i < nf; i++) {
1423 : bates 104 int j, ki = Gp[i+1] - Gp[i],
1424 : bates 105 nci = nc[i], ncip1 = nci + 1, ncisq = nci * nci,
1425 : bates 100 mi = ki/nci;
1426 :     double
1427 : bates 105 *chol = Memcpy(Calloc(ncisq, double),
1428 :     REAL(VECTOR_ELT(Omega, i)), ncisq),
1429 :     *tmp = Calloc(ncisq, double);
1430 :    
1431 : bates 100
1432 : bates 105 F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1433 : bates 104 if (j)
1434 :     error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1435 : bates 105 Memcpy(tmp, chol, ncisq);
1436 : bates 104 F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1437 :     if (j)
1438 :     error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
1439 : bates 100 alpha = (double) -mi;
1440 :     F77_CALL(dsyrk)("U", "N", &nci, &ki,
1441 :     &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1442 : bates 104 &alpha, tmp, &nci);
1443 : bates 176 alpha = (double)(REML ? (nobs-p) : nobs);
1444 : bates 100 F77_CALL(dsyrk)("U", "N", &nci, &mi,
1445 :     &alpha, b + Gp[i], &nci,
1446 : bates 104 &one, tmp, &nci);
1447 : bates 100 if (REML) {
1448 :     for (j = 0; j < p; j++) {
1449 :     F77_CALL(dsyrk)("U", "N", &nci, &mi,
1450 :     &one, RZX + Gp[i] + j*n, &nci,
1451 : bates 104 &one, tmp, &nci);
1452 : bates 100 }
1453 :     }
1454 : bates 104 if (nci == 1) {
1455 :     REAL(ans)[cind++] = *tmp *
1456 :     (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1457 :     } else {
1458 : bates 106 int k, odind = cind + nci;
1459 : bates 104 if (uncst) {
1460 : bates 106 int ione = 1, kk;
1461 : bates 176 double *rr = Calloc(nci, double); /* j'th row of R, the Cholesky factor */
1462 : bates 105 nlme_symmetrize(tmp, nci);
1463 :     for (j = 0; j < nci; j++, cind++) {
1464 : bates 176 for (k = 0; k < j; k++) rr[k] = 0.;
1465 :     for (k = j; k < nci; k++) rr[k] = chol[j + k*nci];
1466 : bates 106 REAL(ans)[cind] = 0.;
1467 : bates 105 for (k = j; k < nci; k++) {
1468 :     for (kk = j; kk < nci; kk++) {
1469 : bates 106 REAL(ans)[cind] += rr[k] * rr[kk] *
1470 : bates 105 tmp[kk * nci + k];
1471 :     }
1472 :     }
1473 :     for (k = j + 1; k < nci; k++) {
1474 : bates 176 REAL(ans)[odind++] = 2. * rr[j] *
1475 :     F77_CALL(ddot)(&nci, rr, &ione, tmp + k*nci, &ione);
1476 : bates 105 }
1477 :     }
1478 : bates 106 Free(rr);
1479 : bates 101 } else {
1480 : bates 104 for (j = 0; j < nci; j++) {
1481 :     REAL(ans)[cind++] = tmp[j * ncip1];
1482 :     for (k = j + 1; k < nci; k++) {
1483 :     REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1484 :     }
1485 :     }
1486 : bates 101 }
1487 : bates 105 cind = odind;
1488 : bates 101 }
1489 : bates 105 Free(tmp); Free(chol);
1490 : bates 100 }
1491 :     UNPROTECT(1);
1492 :     return ans;
1493 :     }
1494 :    
1495 : bates 176 /**
1496 : bates 237 * Return the Hessian of the ML or REML deviance. This is a
1497 :     * placeholder until I work out the evaluation of the analytic
1498 :     * Hessian, which probably will involve several helper functions.
1499 :     *
1500 :     * @param x pointer to an ssclme object
1501 :     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1502 :     * @param Uncp pointer to a logical scalar indicating if the
1503 :     * unconstrained parameterization is to be used
1504 :     *
1505 :     * @return pointer to an approximate Hessian matrix
1506 :     */
1507 :     SEXP ssclme_Hessian(SEXP x, SEXP REMLp, SEXP Uncp)
1508 :     {
1509 :     int j, ncoef = coef_length(length(GET_SLOT(x, Matrix_OmegaSym)),
1510 :     INTEGER(GET_SLOT(x, Matrix_ncSym))),
1511 :     unc = asLogical(Uncp);
1512 :     SEXP ans = PROTECT(allocMatrix(REALSXP, ncoef, ncoef)),
1513 : bates 247 base = PROTECT(ssclme_coef(x, Uncp)),
1514 : bates 237 current = PROTECT(duplicate(base)),
1515 :     gradient;
1516 :    
1517 :     for (j = 0; j < ncoef; j++) {
1518 :     double delta = (REAL(base)[j] ? 1.e-7 * REAL(base)[j] : 1.e-7);
1519 :     int i;
1520 :    
1521 :     for (i = 0; i < ncoef; i++) REAL(current)[i] = REAL(base)[i];
1522 :     REAL(current)[j] += delta/2.;
1523 : bates 247 ssclme_coefGets(x, current, Uncp);
1524 : bates 237 PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1525 :     for (i = 0; i < ncoef; i++) REAL(ans)[j * ncoef + i] = REAL(gradient)[i];
1526 :     UNPROTECT(1);
1527 :     REAL(current)[j] -= delta;
1528 : bates 247 ssclme_coefGets(x, current, Uncp);
1529 : bates 237 PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1530 :     for (i = 0; i < ncoef; i++)
1531 : bates 247 REAL(ans)[j * ncoef + i] =
1532 :     (REAL(ans)[j * ncoef + i] - REAL(gradient)[i])/ delta;
1533 : bates 237 UNPROTECT(1);
1534 : bates 247 for (i = 0; i < j; i++) { /* symmetrize */
1535 : bates 238 REAL(ans)[j * ncoef + i] = REAL(ans)[i * ncoef + j] =
1536 :     (REAL(ans)[j * ncoef + i] + REAL(ans)[i * ncoef + j])/2.;
1537 :     }
1538 : bates 237 }
1539 :     UNPROTECT(3);
1540 :     return ans;
1541 :     }
1542 :    
1543 :     /**
1544 : bates 176 * Calculate and return the fitted values.
1545 :     *
1546 :     * @param x pointer to an ssclme object
1547 :     * @param facs list of grouping factors
1548 :     * @param mmats list of model matrices
1549 :     * @param useRf pointer to a logical scalar indicating if the random effects should be used
1550 :     *
1551 :     * @return pointer to a numeric array of fitted values
1552 :     */
1553 : bates 164 SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats, SEXP useRf)
1554 : bates 28 {
1555 :     SEXP val, b;
1556 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1557 :     i, ione = 1, nf = length(facs), nobs, p;
1558 :     double *vv, one = 1.0, zero = 0.0;
1559 :    
1560 :     if (nf < 1)
1561 :     error("null factor list passed to ssclme_fitted");
1562 :     nobs = length(VECTOR_ELT(facs, 0));
1563 :     val = PROTECT(allocVector(REALSXP, nobs));
1564 :     vv = REAL(val);
1565 :     p = nc[nf] - 1;
1566 :     if (p > 0) {
1567 :     F77_CALL(dgemm)("N", "N", &nobs, &ione, &p, &one,
1568 :     REAL(VECTOR_ELT(mmats, nf)), &nobs,
1569 :     REAL(PROTECT(ssclme_fixef(x))), &p,
1570 :     &zero, vv, &nobs);
1571 :     UNPROTECT(1);
1572 :     } else {
1573 :     memset(vv, 0, sizeof(double) * nobs);
1574 :     }
1575 : bates 164 if (asLogical(useRf)) {
1576 :     b = PROTECT(ssclme_ranef(x));
1577 :     for (i = 0; i < nf; i++) {
1578 :     int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];
1579 :     double *bb = REAL(VECTOR_ELT(b, i)),
1580 :     *mm = REAL(VECTOR_ELT(mmats, i));
1581 :     for (j = 0; j < nobs; ) {
1582 :     int nn = 1, lev = ff[j];
1583 :     /* check for adjacent rows with same factor level */
1584 : bates 282 while ((j + nn) < nobs && ff[j + nn] == lev) nn++;
1585 : bates 164 F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
1586 :     &one, mm + j, &nobs,
1587 :     bb + (lev - 1) * nci, &nci,
1588 :     &one, vv + j, &nobs);
1589 :     j += nn;
1590 :     }
1591 : bates 28 }
1592 : bates 164 UNPROTECT(1);
1593 : bates 28 }
1594 : bates 164 UNPROTECT(1);
1595 : bates 28 return val;
1596 :     }
1597 : bates 45
1598 : bates 122 /**
1599 :     * Return the unscaled variances
1600 :     *
1601 :     * @param x pointer to an ssclme object
1602 :     *
1603 :     * @return a list similar to the Omega list with the unscaled variances
1604 :     */
1605 :     SEXP ssclme_variances(SEXP x)
1606 : bates 45 {
1607 : bates 122 SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1608 : bates 45 int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1609 : bates 57 i, nf = length(Omg);
1610 : bates 45
1611 :     for (i = 0; i < nf; i++) {
1612 : bates 57 double *mm = REAL(VECTOR_ELT(Omg, i));
1613 : bates 122 int j, nci = nc[i];
1614 : bates 45
1615 :     F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1616 :     if (j) /* shouldn't happen */
1617 :     error("DPOTRF returned error code %d on Omega[%d]",
1618 :     j, i + 1);
1619 :     F77_CALL(dpotri)("U", &nci, mm, &nci, &j);
1620 :     if (j) /* shouldn't happen */
1621 :     error("DTRTRI returned error code %d on Omega[%d]",
1622 :     j, i + 1);
1623 : bates 122 nlme_symmetrize(mm, nci);
1624 : bates 45 }
1625 : bates 122 UNPROTECT(1);
1626 :     return Omg;
1627 : bates 45 }
1628 : bates 100
1629 : bates 176 /**
1630 :     * Copy an ssclme object collapsing the fixed effects slots to the response only.
1631 :     *
1632 :     * @param x pointer to an ssclme object
1633 :     *
1634 :     * @return a duplicate of x with the fixed effects slots collapsed to the response only
1635 :     */
1636 : bates 114 SEXP ssclme_collapse(SEXP x)
1637 :     {
1638 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),
1639 :     Omega = GET_SLOT(x, Matrix_OmegaSym),
1640 :     Dim = GET_SLOT(x, Matrix_DimSym);
1641 : bates 164 int nf = length(Omega), nz = INTEGER(Dim)[1];
1642 : bates 114
1643 : bates 176 slot_dup(ans, x, Matrix_DSym);
1644 :     slot_dup(ans, x, Matrix_DIsqrtSym);
1645 :     slot_dup(ans, x, Matrix_DimSym);
1646 :     slot_dup(ans, x, Matrix_GpSym);
1647 :     slot_dup(ans, x, Matrix_LiSym);
1648 :     slot_dup(ans, x, Matrix_LpSym);
1649 :     slot_dup(ans, x, Matrix_LxSym);
1650 :     slot_dup(ans, x, Matrix_OmegaSym);
1651 :     slot_dup(ans, x, Matrix_ParentSym);
1652 :     slot_dup(ans, x, Matrix_bVarSym);
1653 :     slot_dup(ans, x, Matrix_devianceSym);
1654 :     slot_dup(ans, x, Matrix_devCompSym);
1655 :     slot_dup(ans, x, Matrix_iSym);
1656 :     slot_dup(ans, x, Matrix_ncSym);
1657 :     slot_dup(ans, x, Matrix_statusSym);
1658 :     slot_dup(ans, x, Matrix_pSym);
1659 :     slot_dup(ans, x, Matrix_xSym);
1660 : bates 114 INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1661 :     SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));
1662 : bates 122 REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;
1663 : bates 114 SET_SLOT(ans, Matrix_RXXSym, allocMatrix(REALSXP, 1, 1));
1664 : bates 122 REAL(GET_SLOT(ans, Matrix_RXXSym))[0] = NA_REAL;
1665 : bates 114 SET_SLOT(ans, Matrix_ZtXSym, allocMatrix(REALSXP, nz, 1));
1666 :     SET_SLOT(ans, Matrix_RZXSym, allocMatrix(REALSXP, nz, 1));
1667 : bates 115 LOGICAL(GET_SLOT(ans, Matrix_statusSym))[0] = 0;
1668 : bates 114 UNPROTECT(1);
1669 :     return ans;
1670 :     }
1671 :    
1672 : bates 176
1673 :     /**
1674 :     * Create an lme object from its components. This is not done by
1675 :     * new("lme", ...) at the R level because of the possibility of
1676 :     * causing the copying of very large objects.
1677 :     *
1678 :     * @param call Pointer to the original call
1679 :     * @param facs pointer to the list of grouping factors
1680 :     * @param x pointer to the model matrices (may be of length zero)
1681 :     * @param model pointer to the model frame
1682 :     * @param REML pointer to a logical scalar indicating if REML is used
1683 :     * @param rep pointer to the converged ssclme object
1684 :     * @param fitted pointer to the fitted values
1685 :     * @param residuals pointer to the residuals
1686 :     *
1687 :     * @return an lme object
1688 :     */
1689 : bates 164 SEXP ssclme_to_lme(SEXP call, SEXP facs, SEXP x, SEXP model, SEXP REML,
1690 :     SEXP rep, SEXP fitted, SEXP residuals)
1691 :     {
1692 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));
1693 :    
1694 :     SET_SLOT(ans, install("call"), call);
1695 :     SET_SLOT(ans, install("facs"), facs);
1696 :     SET_SLOT(ans, Matrix_xSym, x);
1697 :     SET_SLOT(ans, install("model"), model);
1698 :     SET_SLOT(ans, install("REML"), REML);
1699 :     SET_SLOT(ans, install("rep"), rep);
1700 :     SET_SLOT(ans, install("fitted"), fitted);
1701 :     SET_SLOT(ans, install("residuals"), residuals);
1702 :     UNPROTECT(1);
1703 :     return ans;
1704 :     }
1705 : bates 245

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