# SCM Repository

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

# Annotation of /pkg/src/ssclme.c

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

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