# SCM Repository

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

# Annotation of /pkg/src/ssclme.c

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

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

 root@r-forge.r-project.org ViewVC Help Powered by ViewVC 1.0.0
Thanks to: