SCM

SCM Repository

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

Annotation of /pkg/src/bCrosstab.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge