SCM

SCM Repository

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

Annotation of /pkg/src/sscCrosstab.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 10 #include "sscCrosstab.h"
2 :    
3 :     SEXP sscCrosstab(SEXP flist, SEXP upper)
4 :     {
5 : bates 146 int **fpt, *Ap, *Gp, *TTi, *Ti, *Tj, *dims, i,
6 : bates 10 ncol = 0,
7 :     nfac = length(flist),
8 :     nfc2 = (nfac * (nfac - 1))/2, /* nfac choose 2 */
9 :     nobs = length(VECTOR_ELT(flist, 0)),
10 : bates 146 ntrpl, nz, pos, up = asLogical(upper);
11 :     double *TTx, *Tx;
12 : bates 10 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscCrosstab")));
13 :    
14 :     if (!isNewList(flist) || nfac < 1)
15 : bates 582 error(_("flist must be a non-empty list"));
16 : bates 10 SET_SLOT(val, Matrix_GpSym, allocVector(INTSXP, nfac + 1));
17 :     Gp = INTEGER(GET_SLOT(val, Matrix_GpSym));
18 :     fpt = (int **) R_alloc(nfac, sizeof(int *));
19 :     for (i = 0; i < nfac; i++) {
20 :     SEXP el = VECTOR_ELT(flist, i);
21 :     if (!inherits(el, "factor"))
22 : bates 582 error(_("flist must be a non-empty list of factors"));
23 : bates 10 if (length(el) != nobs)
24 : bates 582 error(_("All elements of flist must have the same length"));
25 : bates 10 Gp[i] = ncol;
26 :     ncol += length(getAttrib(el, R_LevelsSymbol));
27 :     fpt[i] = INTEGER(el);
28 :     }
29 :     Gp[nfac] = ncol;
30 :     SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar(up ? "U" : "L")));
31 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
32 :     dims = INTEGER(GET_SLOT(val, Matrix_DimSym));
33 :     dims[0] = dims[1] = ncol;
34 :     ntrpl = nfc2 * nobs + ncol;
35 :     Ti = Calloc(ntrpl, int); Tj = Calloc(ntrpl, int); TTi = Calloc(ntrpl, int);
36 :     Tx = Calloc(ntrpl, double); TTx = Calloc(ntrpl, double);
37 :     /* Generate the triplet form of the result */
38 :     for (i = 0; i < ncol; i++) {
39 :     Ti[i] = Tj[i] = i; /* The diagonals - these will store counts */
40 :     Tx[i] = 0.0;
41 :     }
42 :     pos = ncol;
43 :     for (i = 0; i < nobs ; i++) {
44 :     int j, jcol, k;
45 :     for (j = 0; j < nfac; j++) {
46 :     jcol = Gp[j] + fpt[j][i] - 1;
47 :     Tx[jcol] += 1.; /* increment diagonal count */
48 :     for (k = j + 1; k < nfac; k++) { /* off-diagonals */
49 :     int irow = Gp[k] + fpt[k][i] - 1;
50 :     if (up) {
51 :     Ti[pos] = jcol; Tj[pos] = irow;
52 :     } else {
53 :     Tj[pos] = jcol; Ti[pos] = irow;
54 :     }
55 :     Tx[pos] = 1.;
56 :     pos++;
57 :     }
58 :     }
59 :     }
60 :     SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, ncol + 1));
61 :     Ap = INTEGER(GET_SLOT(val, Matrix_pSym));
62 : bates 488 triplet_to_col(ncol, ncol, ntrpl, Ti, Tj, Tx, Ap, TTi, TTx);
63 : bates 10 nz = Ap[ncol]; /* non-zeros in Z'Z crosstab */
64 :     SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nz));
65 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nz));
66 :     Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), TTi, nz);
67 :     Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), TTx, nz);
68 :     Free(Ti); Free(Tj); Free(Tx);
69 :     Free(TTi); Free(TTx);
70 : maechler 534
71 : bates 10 UNPROTECT(1);
72 :     return val;
73 :     }
74 : bates 73
75 : bates 143 static
76 : bates 148 void col_metis_order(int j0, int j1, int i2,
77 :     const int Tp[], const int Ti[], int ans[])
78 : bates 143 {
79 : bates 148 int j, nz = 0; /* count off-diagonal pairs */
80 :     for (j = j0; j < j1; j++) { /* columns of interest */
81 :     int ii, nr = 0, p2 = Tp[j + 1];
82 :     for (ii = Tp[j]; ii < p2; ii++) {
83 :     int i = Ti[ii];
84 :     if (j1 <= i && i < i2) nr++; /* verify row index */
85 : bates 143 }
86 : bates 148 nz += (nr * (nr - 1))/2; /* add number of pairs of rows */
87 : bates 143 }
88 : bates 148 if (nz > 0) { /* Form an ssc Matrix */
89 :     int j, n = i2 - j1, /* number of rows */
90 :     nnz = n + nz, pos;
91 :     int *Ap = Calloc(n + 1, int),
92 :     *Ai = Calloc(nnz, int),
93 :     *Tj = Calloc(nnz, int),
94 : bates 149 *TTi = Calloc(nnz, int),
95 :     *perm = Calloc(n, int),
96 :     *iperm = Calloc(n, int);
97 : bates 143
98 : bates 148 for (j = 0; j < n; j++) { /* diagonals */
99 :     TTi[j] = Tj[j] = j;
100 : bates 85 }
101 : bates 148 pos = n;
102 :     for (j = j0; j < j1; j++) { /* create the pairs */
103 : bates 153 int ii, p2 = Tp[j + 1];
104 : bates 148 for (ii = Tp[j]; ii < p2; ii++) {
105 :     int r1 = Ti[ii], i1;
106 :     if (j1 <= r1 && r1 < i2) {
107 :     for (i1 = ii + 1; i1 < p2; i1++) {
108 :     int r2 = Ti[i1];
109 :     if (r2 < i2) {
110 :     TTi[pos] = r2 - j1;
111 :     Tj[pos] = r1 - j1;
112 :     pos++;
113 :     }
114 :     }
115 : bates 93 }
116 : bates 85 }
117 :     }
118 : bates 488 triplet_to_col(n, n, nnz, TTi, Tj, (double *) NULL,
119 : bates 153 Ap, Ai, (double *) NULL);
120 : bates 148 ssc_metis_order(n, Ap, Ai, perm, iperm);
121 :     for (j = j1; j < i2; j++) ans[j] = j1 + iperm[j - j1];
122 : bates 149 Free(TTi); Free(Tj); Free(Ai); Free(Ap);
123 : bates 148 Free(perm); Free(iperm);
124 : bates 85 }
125 :     }
126 :    
127 : bates 127 SEXP sscCrosstab_groupedPerm(SEXP ctab)
128 : bates 74 {
129 :     SEXP
130 :     GpSlot = GET_SLOT(ctab, Matrix_GpSym),
131 :     iSlot = GET_SLOT(ctab, Matrix_iSym),
132 :     pSlot = GET_SLOT(ctab, Matrix_pSym);
133 :     int *Ai = INTEGER(iSlot),
134 :     *Ap = INTEGER(pSlot),
135 :     *Gp = INTEGER(GpSlot),
136 : bates 85 i,
137 : bates 74 n = length(pSlot) - 1, /* number of columns */
138 :     nf = length(GpSlot) - 1, /* number of factors */
139 : bates 143 up;
140 : bates 74 SEXP ans = PROTECT(allocVector(INTSXP, n));
141 :    
142 : maechler 534 up = *CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0)) != 'L';
143 : bates 143 if (nf > 1 && up) { /* transpose */
144 :     int nz = length(iSlot);
145 :     int *ai = Calloc(nz, int),
146 :     *ap = Calloc(n + 1, int);
147 :     double *ax = Calloc(nz, double);
148 : bates 74
149 : bates 143 csc_components_transpose(n, n, nz, Ap, Ai,
150 :     REAL(GET_SLOT(ctab, Matrix_xSym)),
151 :     ap, ai, ax);
152 :     Ap = ap;
153 :     Ai = ai;
154 :     Free(ax); /* don't need values, only positions */
155 :     }
156 : bates 85 for (i = 0; i < n; i++) {
157 : bates 86 INTEGER(ans)[i] = i; /* initialize permutation to identity */
158 : bates 85 }
159 : bates 146 for (i = 1; i < nf; i++) {
160 :     col_metis_order(Gp[i - 1], Gp[i], Gp[i+1], Ap, Ai, INTEGER(ans));
161 :     }
162 : bates 143 if (nf > 1 && up) {Free(Ap); Free(Ai);}
163 : bates 74 UNPROTECT(1);
164 :     return ans;
165 :     }
166 : bates 153
167 : maechler 534 /**
168 : bates 153 * Project the (2,1) component of an sscCrosstab object into the (2,2)
169 :     * component (for illustration only)
170 : maechler 534 *
171 : bates 153 * @param ctab pointer to a sscCrosstab object
172 : maechler 534 *
173 : bates 478 * @return a pointer to an dsCMatrix giving the projection of the 2,1 component
174 : bates 153 */
175 :     SEXP sscCrosstab_project(SEXP ctab)
176 :     {
177 :     SEXP
178 :     GpSlot = GET_SLOT(ctab, Matrix_GpSym),
179 :     iSlot = GET_SLOT(ctab, Matrix_iSym),
180 :     pSlot = GET_SLOT(ctab, Matrix_pSym);
181 :     int *Ai = INTEGER(iSlot),
182 :     *Ap = INTEGER(pSlot),
183 :     *Gp = INTEGER(GpSlot),
184 : bates 174 j, j0, j1, i2,
185 : bates 153 n = length(pSlot) - 1, /* number of columns */
186 :     nf = length(GpSlot) - 1, /* number of factors */
187 :     nz, up;
188 : bates 478 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix")));
189 : bates 153
190 : maechler 534 up = *CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0)) != 'L';
191 : bates 153 if (nf > 1 && up) { /* transpose */
192 :     int nz = length(iSlot);
193 :     int *ai = Calloc(nz, int),
194 :     *ap = Calloc(n + 1, int);
195 :     double *ax = Calloc(nz, double);
196 :    
197 :     csc_components_transpose(n, n, nz, Ap, Ai,
198 :     REAL(GET_SLOT(ctab, Matrix_xSym)),
199 :     ap, ai, ax);
200 :     Ap = ap;
201 :     Ai = ai;
202 :     Free(ax); /* don't need values, only positions */
203 :     }
204 : maechler 534
205 : bates 153 nz = 0; /* count of off-diagonal pairs */
206 :     j0 = 0; j1 = Gp[1]; i2 = Gp[2];
207 :     for (j = j0; j < j1; j++) { /* columns of interest */
208 :     int ii, nr = 0, p2 = Ap[j + 1];
209 :     for (ii = Ap[j]; ii < p2; ii++) {
210 :     int i = Ai[ii];
211 :     if (j1 <= i && i < i2) nr++; /* verify row index */
212 :     }
213 :     nz += (nr * (nr - 1))/2; /* add number of pairs of rows */
214 :     }
215 :     if (nz > 0) { /* Form an ssc Matrix */
216 :     int j, n = i2 - j1, /* number of rows */
217 :     nnz = n + nz, pos;
218 :     int *AAp,
219 :     *AAi = Calloc(nnz, int),
220 :     *Tj = Calloc(nnz, int),
221 :     *TTi = Calloc(nnz, int);
222 :     double *Ax;
223 :    
224 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));
225 :     AAp = INTEGER(GET_SLOT(ans, Matrix_pSym));
226 :     for (j = 0; j < n; j++) { /* diagonals */
227 :     TTi[j] = Tj[j] = j;
228 :     }
229 :     pos = n;
230 :     for (j = j0; j < j1; j++) { /* create the pairs */
231 :     int ii, p2 = Ap[j + 1];
232 :     for (ii = Ap[j]; ii < p2; ii++) {
233 :     int r1 = Ai[ii], i1;
234 :     if (j1 <= r1 && r1 < i2) {
235 :     for (i1 = ii + 1; i1 < p2; i1++) {
236 :     int r2 = Ai[i1];
237 :     if (r2 < i2) {
238 :     TTi[pos] = r2 - j1;
239 :     Tj[pos] = r1 - j1;
240 :     pos++;
241 :     }
242 :     }
243 :     }
244 :     }
245 :     }
246 : bates 488 triplet_to_col(n, n, nnz, TTi, Tj, (double *) NULL,
247 : bates 153 AAp, AAi, (double *) NULL);
248 :     nz = AAp[n];
249 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nz));
250 :     Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), AAi, nz);
251 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nz));
252 :     Ax = REAL(GET_SLOT(ans, Matrix_xSym));
253 :     for (j = 0; j < nz; j++) Ax[j] = 1.;
254 : maechler 534 SET_SLOT(ans, Matrix_uploSym, mkString("L"));
255 : bates 153 SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
256 :     AAp = INTEGER(GET_SLOT(ans, Matrix_DimSym));
257 :     AAp[0] = AAp[1] = n;
258 :     Free(TTi); Free(Tj); Free(AAi);
259 :     }
260 :     if (nf > 1 && up) {Free(Ap); Free(Ai);}
261 :     UNPROTECT(1);
262 :     return ans;
263 :     }
264 : bates 184
265 : maechler 534 /**
266 : bates 184 * Project the first group of columns in an sscCrosstab object onto the
267 :     * remaining columns.
268 : maechler 534 *
269 : bates 184 * @param ctab pointer to a sscCrosstab object
270 : maechler 534 *
271 : bates 478 * @return a pointer to an dsCMatrix with the projection
272 : bates 184 */
273 :     SEXP sscCrosstab_project2(SEXP ctab)
274 :     {
275 :     SEXP
276 :     GpSlot = GET_SLOT(ctab, Matrix_GpSym),
277 :     iSlot = GET_SLOT(ctab, Matrix_iSym),
278 :     pSlot = GET_SLOT(ctab, Matrix_pSym);
279 :     int *Ai = INTEGER(iSlot),
280 :     *Ap = INTEGER(pSlot),
281 :     *Gp = INTEGER(GpSlot),
282 :     i, i2, j, j1, k, k2,
283 :     nf = length(GpSlot) - 1, /* number of factors */
284 :     nz, pos, up, *AAp, *Ti, *Cp, *ind;
285 :     double *Ax = REAL(GET_SLOT(ctab, Matrix_xSym));
286 : bates 478 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix")));
287 : bates 184
288 : maechler 534
289 : bates 582 if (nf < 2) error(_("sscCrosstab_project2 requires more than one group"));
290 : maechler 534 up = *CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0)) != 'L';
291 : bates 184 if (up) { /* tranpose */
292 :     int n = length(pSlot) - 1; /* number of columns */
293 :     int nnz = length(iSlot);
294 :     int *ai = Calloc(nnz, int);
295 :     int *ap = Calloc(n + 1, int);
296 :     double *ax = Calloc(nnz, double);
297 :    
298 :     csc_components_transpose(n, n, nnz, Ap, Ai, Ax, ap, ai, ax);
299 :     Ap = ap; Ai = ai; Ax = ax;
300 :     }
301 :    
302 :     j1 = Gp[1]; i2 = Gp[nf]; /* group boundaries */
303 :     Cp = Calloc(j1, int); /* first pos in col with row > i */
304 :     for (j = 0; j < j1; j++) Cp[j] = Ap[j] + 1; /* Ap[j] is the diagonal */
305 :    
306 :     nz = Ap[i2] - Ap[j1]; /* upper bound on nonzeros in result */
307 : maechler 534 for (j = 0; j < j1; j++) {
308 : bates 184 int nr = Ap[j + 1] - Ap[j] - 1;
309 :     nz += (nr * (nr - 1))/2; /* add number of pairs of rows below diag */
310 :     }
311 :    
312 :     Ti = Calloc(nz, int); /* temporary row indices */
313 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, i2 - j1 + 1));
314 :     AAp = INTEGER(GET_SLOT(ans, Matrix_pSym)); /* column pointers */
315 :    
316 :     AAp[0] = 0; pos = 0;
317 :     ind = Calloc(i2 - j1, int); /* indicator of rows in same column */
318 :     for (i = j1; i < i2; i++) {
319 :     for (k = j1; k < i2; k++) ind[k - j1] = 0;
320 :     for (j = 0; j < j1; j++) {
321 :     if (Ai[Cp[j]] == i) { /* go down the column */
322 :     k2 = Ap[j+1];
323 :     for(k = Cp[j] + 1; k < k2; k++) ind[Ai[k] - j1] = 1;
324 :     Cp[j]++;
325 :     }
326 :     }
327 : maechler 534
328 : bates 184 Ti[pos++] = i - j1; /* diagonal element */
329 :     for (k = i+1; k < i2; k++) { /* projected pairs */
330 :     int ii = k - j1;
331 :     if (ind[ii]) Ti[pos++] = ii;
332 :     }
333 :     k2 = Ap[i+1];
334 :     for (k = Ap[i] + 1; k < k2; k++) { /* previous off-diagonals */
335 :     Ti[pos++] = Ai[k] - j1;
336 :     }
337 :     AAp[i - j1 + 1] = pos;
338 :     }
339 :     nz = AAp[i2 - j1];
340 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nz));
341 :     Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), Ti, nz);
342 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nz));
343 :     Ax = REAL(GET_SLOT(ans, Matrix_xSym));
344 :     for (j = 0; j < nz; j++) Ax[j] = 1.;
345 : maechler 534 SET_SLOT(ans, Matrix_uploSym, mkString("L"));
346 : bates 184 SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
347 :     AAp = INTEGER(GET_SLOT(ans, Matrix_DimSym));
348 :     AAp[0] = AAp[1] = i2 - j1;
349 :     Free(Ti); Free(Cp); Free(ind);
350 :     if (up) {Free(Ap); Free(Ai); free(Ax);}
351 :     UNPROTECT(1);
352 :     return ans;
353 :     }

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