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 159 - (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 :     error("flist must be a non-empty list");
16 :     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 :     error("flist must be a non-empty list of factors");
23 :     if (length(el) != nobs)
24 :     error("All elements of flist must have the same length");
25 :     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 :     triplet_to_col(ncol, ncol, ntrpl, Ti, Tj, Tx, Ap, TTi, TTx);
63 :     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 :    
71 :     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 153 triplet_to_col(n, n, nnz, TTi, Tj, (double *) NULL,
119 :     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 : bates 143 up = toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L';
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 :     /**
168 :     * Project the (2,1) component of an sscCrosstab object into the (2,2)
169 :     * component (for illustration only)
170 :     *
171 :     * @param ctab pointer to a sscCrosstab object
172 :     *
173 :     * @return a pointer to an tscMatrix giving the projection of the 2,1 component
174 :     */
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 :     i, j, j0, j1, i2,
185 :     n = length(pSlot) - 1, /* number of columns */
186 :     nf = length(GpSlot) - 1, /* number of factors */
187 :     nz, up;
188 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix")));
189 :    
190 :     up = toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L';
191 :     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 :    
205 :     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 :     triplet_to_col(n, n, nnz, TTi, Tj, (double *) NULL,
247 :     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 :     SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));
255 :     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 :     }

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