# SCM Repository

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

# Annotation of /pkg/src/bCrosstab.c

 1 : bates 356 /* TODO: 1) combine diag_update and offdiag_update - DONE 2 : 2) Move Parent_inverse from tscMatrix to here 3 : 3) Modify the creation of the L diagonals to include a call 4 : to ldl_numeric to get Li. 5 : */ 6 : 7 : bates 347 #include "bCrosstab.h" 8 : 9 : /** 10 : * Calculate the zero-based index in a row-wise packed lower 11 : * triangular matrix. This is used for the arrays of blocked sparse matrices. 12 : * 13 : * @param i row number (0-based) 14 : * @param k column number (0-based) 15 : * 16 : * @return The 0-based index of the (i,k) element of a row-wise packed lower 17 : * triangular matrix. 18 : */ 19 : static R_INLINE 20 : int Lind(int i, int k) 21 : { 22 : return (i * (i + 1))/2 + k; 23 : } 24 : 25 : /** 26 : * Permute an index vector 27 : * 28 : * @param i vector of 0-based indices 29 : * @param nnz length of vector i 30 : * @param perm 0-based permutation vector of length max(i) 31 : */ 32 : static R_INLINE 33 : bates 356 void ind_permute(int i[], int nnz, const int perm[]) 34 : bates 347 { 35 : int j; 36 : for (j = 0; j < nnz; j++) i[j] = perm[i[j]]; 37 : } 38 : 39 : /** 40 : bates 356 * Force indices to be in the upper triangle of a matrix 41 : bates 347 * 42 : * @param i vector of 0-based row indices 43 : * @param j vector of 0-based column indices 44 : * @param nnz length of index vectors 45 : */ 46 : static R_INLINE 47 : bates 356 void make_upper_triangular(int i[], int j[], int nnz) 48 : bates 347 { 49 : int k; 50 : for (k = 0; k < nnz; k++) { 51 : bates 356 if (i[k] > j[k]) { 52 : bates 347 int tmp = i[k]; 53 : i[k] = j[k]; 54 : j[k] = tmp; 55 : } 56 : } 57 : } 58 : 59 : /** 60 : * Create a named list of length n 61 : * 62 : * @param n length of list to return 63 : * @param names names of list elements 64 : * 65 : * @return pointer to a named list of length n 66 : */ 67 : static 68 : SEXP make_named_list(int n, char **names) 69 : { 70 : SEXP ans = PROTECT(allocVector(VECSXP, n)), 71 : bates 356 nms = PROTECT(allocVector(STRSXP, n)); 72 : bates 347 int i; 73 : 74 : bates 356 for (i = 0; i < n; i++) SET_STRING_ELT(nms, i, mkChar(names[i])); 75 : bates 347 setAttrib(ans, R_NamesSymbol, nms); 76 : UNPROTECT(2); 77 : return ans; 78 : } 79 : 80 : bates 356 static R_INLINE 81 : int check_csc_index(const int p[], const int i[], int row, int col) 82 : { 83 : int k, k2 = p[col + 1]; 84 : /* use a linear search for now */ 85 : /* perhaps replace by bsearch later */ 86 : for (k = p[col]; k < k2; k++) { 87 : if (i[k] == row) return k; 88 : } 89 : return -1; 90 : } 91 : 92 : bates 347 /** 93 : * Update a diagonal block 94 : * 95 : * @param ctab pointer to a blocked crosstabulation object 96 : bates 356 * @param j index of updating column 97 : bates 347 * @param i index of diagonal block to be updated 98 : */ 99 : bates 356 static void diag_update(SEXP ctab, int j, int i) 100 : bates 347 { 101 : SEXP db = VECTOR_ELT(ctab, Lind(i, i)), 102 : bates 356 jb = VECTOR_ELT(ctab, Lind(i, j)); 103 : bates 347 SEXP dpp = GET_SLOT(db, Matrix_pSym), 104 : jpp = GET_SLOT(jb, Matrix_pSym); 105 : bates 356 int *di = INTEGER(GET_SLOT(db, Matrix_iSym)), 106 : bates 347 *dp = INTEGER(dpp), 107 : bates 356 *ji = INTEGER(GET_SLOT(jb, Matrix_iSym)), 108 : bates 347 *jp = INTEGER(jpp), 109 : bates 356 dnc = length(dpp) - 1, 110 : jnc = length(jpp) - 1; 111 : int jj, extra; 112 : /* bound the number of extra elements */ 113 : extra = 0; 114 : for (jj = 0; jj < jnc; jj++) { 115 : int k, kk, k2 = jp[jj + 1]; 116 : for (k = jp[jj]; k < k2; k++) { 117 : for (kk = k; kk < k2; kk++) { 118 : if (check_csc_index(dp, di, ji[k], ji[kk]) < 0) extra++; 119 : } 120 : } 121 : } 122 : if (!extra) return; 123 : { 124 : int pos, nnz = dp[dnc]; 125 : int ntot = nnz + extra; 126 : int *Ai = Calloc(ntot, int), 127 : *Ti = Calloc(ntot, int), 128 : *Tj = Calloc(ntot, int); 129 : double *Ax; 130 : 131 : Memcpy(Ti, di, nnz); /* make a copy of the row indices */ 132 : pos = 0; /* fill in the column indices */ 133 : for (jj = 0; jj < dnc; jj++) { 134 : int j2 = dp[jj + 1]; 135 : while (pos < j2) { 136 : Tj[pos] = jj; 137 : pos++; 138 : } 139 : } 140 : /* add the extra elements */ 141 : for (jj = 0; jj < jnc; jj++) { 142 : int k, kk, k2 = jp[jj + 1]; 143 : for (k = jp[jj]; k < k2; k++) { 144 : for (kk = k; kk < k2; kk++) { 145 : if (check_csc_index(dp, di, ji[k], ji[kk]) < 0) { 146 : Ti[pos] = ji[k]; 147 : Tj[pos] = ji[kk]; 148 : pos++; 149 : } 150 : } 151 : } 152 : } 153 : triplet_to_col(dnc, dnc, ntot, Ti, Tj, (double *) NULL, 154 : dp, Ai, (double *) NULL); 155 : nnz = dp[dnc]; 156 : SET_SLOT(db, Matrix_iSym, allocVector(INTSXP, nnz)); 157 : Memcpy(INTEGER(GET_SLOT(db, Matrix_iSym)), Ai, nnz); 158 : SET_SLOT(db, Matrix_xSym, allocVector(REALSXP, nnz)); 159 : Ax = REAL(GET_SLOT(db, Matrix_xSym)); 160 : for (j = 0; j < nnz; j++) Ax[j] = 1.; 161 : Free(Ai); Free(Ti); Free(Tj); 162 : return; 163 : } 164 : bates 347 } 165 : 166 : /** 167 : bates 356 * Update a block 168 : bates 347 * 169 : * @param ctab pointer to a blocked crosstabulation object 170 : bates 356 * @param j index of updating column 171 : * @param k column index of block to be updated 172 : * @param i row index of block to be updated (j < k <= i) 173 : bates 347 */ 174 : bates 356 static void block_update(SEXP ctab, int j, int k, int i) 175 : bates 347 { 176 : bates 356 SEXP tb = VECTOR_ELT(ctab, Lind(i, k)), 177 : ib = VECTOR_ELT(ctab, Lind(i, j)), 178 : kb = VECTOR_ELT(ctab, Lind(k, j)); 179 : SEXP tpp = GET_SLOT(tb, Matrix_pSym), 180 : kpp = GET_SLOT(kb, Matrix_pSym); 181 : int *ti = INTEGER(GET_SLOT(tb, Matrix_iSym)), 182 : *tp = INTEGER(tpp), 183 : *ii = INTEGER(GET_SLOT(ib, Matrix_iSym)), 184 : *ip = INTEGER(GET_SLOT(ib, Matrix_pSym)), 185 : *ki = INTEGER(GET_SLOT(kb, Matrix_iSym)), 186 : *kp = INTEGER(kpp), 187 : tnc = length(tpp) - 1, 188 : knc = length(kpp) - 1; 189 : int jj, extra; 190 : bates 347 191 : bates 356 if (k > i || j >= k) 192 : error("i,j,k values of %d,%d,%d do not satisfy j < k <= i", 193 : i, j, k); 194 : /* bound the number of extra elements */ 195 : extra = 0; 196 : for (jj = 0; jj < knc; jj++) { 197 : int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1]; 198 : for (kk = kp[jj]; kk < k2; kk++) { 199 : for (i1 = ip[jj]; i1 < i2; i1++) { 200 : if ((check_csc_index(tp, ti, ii[i1], ki[kk]) < 0) && 201 : /* only update upper triangle of 202 : * diagonal blocks */ 203 : ((k != i) || (ii[i1] <= ki[kk]))) extra++; 204 : } 205 : } 206 : bates 347 } 207 : bates 356 if (!extra) return; 208 : { 209 : int pos, nnz = tp[tnc]; 210 : int ntot = nnz + extra; 211 : int *Ai = Calloc(ntot, int), 212 : *Ti = Calloc(ntot, int), 213 : *Tj = Calloc(ntot, int), 214 : *Dims = INTEGER(GET_SLOT(tb, Matrix_DimSym)); 215 : double *Ax; 216 : 217 : Memcpy(Ti, ti, nnz); /* make a copy of the row indices */ 218 : for (pos = 0, jj = 0; jj < tnc; jj++) { /* fill in the column indices */ 219 : int j2 = tp[jj + 1]; 220 : for (; pos < j2; pos++) Tj[pos] = jj; 221 : } 222 : /* add the extra elements */ 223 : for (jj = 0; jj < knc; jj++) { 224 : int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1]; 225 : for (kk = kp[jj]; kk < k2; kk++) { 226 : for (i1 = ip[jj]; i1 < i2; i1++) { 227 : if ((check_csc_index(tp, ti, ii[i1], ki[kk]) < 0) && 228 : ((k != i) || (ii[i1] <= ki[kk]))) { 229 : Ti[pos] = ii[i1]; 230 : Tj[pos] = ki[kk]; 231 : pos++; 232 : } 233 : } 234 : } 235 : } 236 : triplet_to_col(Dims[0], Dims[1], ntot, Ti, Tj, (double *) NULL, 237 : tp, Ai, (double *) NULL); 238 : nnz = tp[tnc]; 239 : SET_SLOT(tb, Matrix_iSym, allocVector(INTSXP, nnz)); 240 : Memcpy(INTEGER(GET_SLOT(tb, Matrix_iSym)), Ai, nnz); 241 : SET_SLOT(tb, Matrix_xSym, allocVector(REALSXP, nnz)); 242 : Ax = REAL(GET_SLOT(tb, Matrix_xSym)); 243 : for (j = 0; j < nnz; j++) Ax[j] = 1.; 244 : Free(Ai); Free(Ti); Free(Tj); 245 : return; 246 : } 247 : bates 347 } 248 : 249 : /** 250 : * Permute the levels of one of the grouping factors in a bCrosstab object 251 : * 252 : * @param ctab Pointer to a bCrosstab object 253 : bates 356 * @param nf number of factors in ctab 254 : * @param jj index (0-based) of the factor levels to permute 255 : * @param ncj number of columns in level jj 256 : bates 347 * @param perm permutation (0-based) to apply 257 : bates 356 * @param pperm inverse of the permutation 258 : bates 347 */ 259 : bates 356 void bCrosstab_permute(SEXP ctab, int nf, int jj, 260 : const int perm[], const int iperm[]) 261 : bates 347 { 262 : bates 356 SEXP trip, ipt, jpt; 263 : int j, nnz; 264 : bates 347 265 : bates 356 for (j = 0; j < nf; j++) { 266 : int ind = (j < jj ? Lind(jj, j) : Lind(j, jj)); 267 : SEXP csc = VECTOR_ELT(ctab, ind); 268 : int *Dims = INTEGER(GET_SLOT(csc, Matrix_DimSym)); 269 : 270 : trip = csc_to_triplet(csc); 271 : bates 347 ipt = GET_SLOT(trip, Matrix_iSym); 272 : bates 356 nnz = length(ipt); 273 : bates 347 jpt = GET_SLOT(trip, Matrix_jSym); 274 : bates 356 if (j <= jj) ind_permute(INTEGER(ipt), nnz, iperm); 275 : if (j >= jj) ind_permute(INTEGER(jpt), nnz, iperm); 276 : if (j == jj) 277 : make_upper_triangular(INTEGER(ipt), INTEGER(jpt), nnz); 278 : triplet_to_col(Dims[0], Dims[1], nnz, INTEGER(ipt), 279 : INTEGER(jpt), 280 : REAL(GET_SLOT(trip, Matrix_xSym)), 281 : INTEGER(GET_SLOT(csc, Matrix_pSym)), 282 : INTEGER(GET_SLOT(csc, Matrix_iSym)), 283 : REAL(GET_SLOT(csc, Matrix_xSym))); 284 : bates 347 } 285 : } 286 : 287 : bates 356 static void 288 : factor_levels_permute(SEXP dest, SEXP src, const int perm[], 289 : const int iperm[]) 290 : { 291 : SEXP dlev = getAttrib(dest, R_LevelsSymbol), 292 : slev = getAttrib(src, R_LevelsSymbol); 293 : int nlev = length(dlev), flen = length(dest); 294 : int *d = INTEGER(dest), *s = INTEGER(src), i; 295 : 296 : if (length(slev) != nlev) 297 : error("number of levels in src and dest must match"); 298 : if (length(src) != flen) 299 : error("length of src and dest must match"); 300 : for (i = 0; i < nlev; i++) 301 : SET_STRING_ELT(dlev, i, STRING_ELT(slev, perm[i])); 302 : for (i = 0; i < flen; i++) 303 : d[i] = 1 + iperm[s[i]-1]; 304 : } 305 : 306 : SEXP bCrosstab_convert(SEXP bCtab) 307 : { 308 : char *anms[] = {"flist", "L", "Linv", "perm", "Parent"}; 309 : SEXP flist = VECTOR_ELT(bCtab, 0), 310 : ctab = PROTECT(duplicate(VECTOR_ELT(bCtab, 1))), 311 : ans = PROTECT(make_named_list(5, anms)); 312 : SEXP fcp, perm, L, Linv, Parent; 313 : int ctbl = length(ctab), j, nf = length(flist); 314 : 315 : if (ctbl != (nf*(nf + 1))/2) 316 : error("length of ctab = %d is not permisable", ctbl); 317 : SET_VECTOR_ELT(ans, 0, duplicate(flist)); 318 : SET_VECTOR_ELT(ans, 1, allocVector(VECSXP, ctbl)); 319 : SET_VECTOR_ELT(ans, 2, allocVector(VECSXP, nf)); 320 : SET_VECTOR_ELT(ans, 3, allocVector(VECSXP, nf)); 321 : SET_VECTOR_ELT(ans, 4, allocVector(VECSXP, nf)); 322 : fcp = VECTOR_ELT(ans, 0); 323 : L = VECTOR_ELT(ans, 1); 324 : setAttrib(L, R_NamesSymbol, duplicate(getAttrib(ctab, R_NamesSymbol))); 325 : Linv = VECTOR_ELT(ans, 2); 326 : setAttrib(Linv, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol))); 327 : perm = VECTOR_ELT(ans, 3); 328 : setAttrib(perm, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol))); 329 : Parent = VECTOR_ELT(ans, 4); 330 : setAttrib(Parent, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol))); 331 : for (j = 0; j < nf; j++) { 332 : int dind = Lind(j, j), i, k; 333 : SEXP ctd = VECTOR_ELT(ctab, dind); /* diagonal in crosstab */ 334 : SEXP Dimslot = GET_SLOT(ctd, Matrix_DimSym), 335 : Ljj, cpp = GET_SLOT(ctd, Matrix_pSym), 336 : cip = GET_SLOT(ctd, Matrix_iSym); 337 : int *Lp, *Perm, *cp = INTEGER(cpp), 338 : *ci = INTEGER(cip), *parent, 339 : ncj = length(cpp) - 1, 340 : nnz = length(cip); 341 : 342 : SET_VECTOR_ELT(Parent, j, allocVector(INTSXP, ncj)); 343 : parent = INTEGER(VECTOR_ELT(Parent, j)); 344 : SET_VECTOR_ELT(perm, j, allocVector(INTSXP, ncj)); 345 : Perm = INTEGER(VECTOR_ELT(perm, j)); 346 : SET_VECTOR_ELT(L, dind, NEW_OBJECT(MAKE_CLASS("cscMatrix"))); 347 : Ljj = VECTOR_ELT(L, dind); 348 : SET_SLOT(Ljj, Matrix_DimSym, duplicate(Dimslot)); 349 : SET_SLOT(Ljj, Matrix_factorization, allocVector(VECSXP, 0)); 350 : SET_SLOT(Ljj, Matrix_pSym, allocVector(INTSXP, ncj + 1)); 351 : Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym)); 352 : if (nnz > ncj) { /* calculate fill-reducing permutation */ 353 : int *iPerm = Calloc(ncj, int), 354 : *Lnz = Calloc(ncj, int); 355 : double *Lx; 356 : 357 : ssc_metis_order(ncj, cp, ci, Perm, iPerm); 358 : /* apply to the crosstabulation */ 359 : bCrosstab_permute(ctab, nf, j, Perm, iPerm); 360 : /* apply to the factor */ 361 : factor_levels_permute(VECTOR_ELT(fcp, j), 362 : VECTOR_ELT(flist, j), Perm, iPerm); 363 : /* symbolic analysis to get Parent */ 364 : ldl_symbolic(ncj, cp, ci, Lp, /* iPerm used as scratch */ 365 : parent, Lnz, iPerm, 366 : (int *) NULL, (int *) NULL); 367 : nnz = Lp[ncj]; 368 : Free(iPerm); Free(Lnz); 369 : SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz)); 370 : SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, nnz)); 371 : Lnz = INTEGER(GET_SLOT(Ljj, Matrix_iSym)); 372 : Lx = REAL(GET_SLOT(Ljj, Matrix_xSym)); 373 : for (i = 0; i < nnz; i++) { 374 : Lnz[i] = 0; Lx[i] = 0.; 375 : } 376 : } else { 377 : for (i = 0; i <= ncj; i++) { 378 : Lp[i] = 0; 379 : parent[i] = -1; 380 : Perm[i] = i; 381 : } 382 : SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0)); 383 : SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, 0)); 384 : } 385 : SET_VECTOR_ELT(Linv, j, 386 : Parent_inverse(VECTOR_ELT(Parent, j), 387 : ScalarLogical(1))); 388 : /* FIXME: Update any blocks below the diagonal in this column if necessary*/ 389 : for (k = j+1; k < nf; k++) { /* update remaining columns, if any */ 390 : for (i = k; i < nf; i++) block_update(ctab, j, k, i); 391 : } 392 : for (i= 0; i < j; i++) { /* copy blocks to the left */ 393 : SET_VECTOR_ELT(L, Lind(j,i), duplicate(VECTOR_ELT(ctab, Lind(j,i)))); 394 : } 395 : } 396 : UNPROTECT(2); 397 : return ans; 398 : } 399 : 400 : 401 : 402 : bates 347