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 86 - (view) (download) (as text)

1 : bates 10 #include "sscCrosstab.h"
2 :    
3 :     SEXP sscCrosstab(SEXP flist, SEXP upper)
4 :     {
5 :     int
6 :     **fpt,
7 :     *Ap,
8 :     *Gp,
9 :     *TTi,
10 :     *Ti,
11 :     *Tj,
12 :     *dims,
13 :     i,
14 :     ncol = 0,
15 :     nfac = length(flist),
16 :     nfc2 = (nfac * (nfac - 1))/2, /* nfac choose 2 */
17 :     nobs = length(VECTOR_ELT(flist, 0)),
18 :     ntrpl,
19 :     nz,
20 :     pos,
21 :     up = asLogical(upper);
22 :     double
23 :     *TTx,
24 :     *Tx;
25 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscCrosstab")));
26 :    
27 :     if (!isNewList(flist) || nfac < 1)
28 :     error("flist must be a non-empty list");
29 :     SET_SLOT(val, Matrix_GpSym, allocVector(INTSXP, nfac + 1));
30 :     Gp = INTEGER(GET_SLOT(val, Matrix_GpSym));
31 :     fpt = (int **) R_alloc(nfac, sizeof(int *));
32 :     for (i = 0; i < nfac; i++) {
33 :     SEXP el = VECTOR_ELT(flist, i);
34 :     if (!inherits(el, "factor"))
35 :     error("flist must be a non-empty list of factors");
36 :     if (length(el) != nobs)
37 :     error("All elements of flist must have the same length");
38 :     Gp[i] = ncol;
39 :     ncol += length(getAttrib(el, R_LevelsSymbol));
40 :     fpt[i] = INTEGER(el);
41 :     }
42 :     Gp[nfac] = ncol;
43 :     SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar(up ? "U" : "L")));
44 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
45 :     dims = INTEGER(GET_SLOT(val, Matrix_DimSym));
46 :     dims[0] = dims[1] = ncol;
47 :     ntrpl = nfc2 * nobs + ncol;
48 :     Ti = Calloc(ntrpl, int); Tj = Calloc(ntrpl, int); TTi = Calloc(ntrpl, int);
49 :     Tx = Calloc(ntrpl, double); TTx = Calloc(ntrpl, double);
50 :     /* Generate the triplet form of the result */
51 :     for (i = 0; i < ncol; i++) {
52 :     Ti[i] = Tj[i] = i; /* The diagonals - these will store counts */
53 :     Tx[i] = 0.0;
54 :     }
55 :     pos = ncol;
56 :     for (i = 0; i < nobs ; i++) {
57 :     int j, jcol, k;
58 :     for (j = 0; j < nfac; j++) {
59 :     jcol = Gp[j] + fpt[j][i] - 1;
60 :     Tx[jcol] += 1.; /* increment diagonal count */
61 :     for (k = j + 1; k < nfac; k++) { /* off-diagonals */
62 :     int irow = Gp[k] + fpt[k][i] - 1;
63 :     if (up) {
64 :     Ti[pos] = jcol; Tj[pos] = irow;
65 :     } else {
66 :     Tj[pos] = jcol; Ti[pos] = irow;
67 :     }
68 :     Tx[pos] = 1.;
69 :     pos++;
70 :     }
71 :     }
72 :     }
73 :     SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, ncol + 1));
74 :     Ap = INTEGER(GET_SLOT(val, Matrix_pSym));
75 :     triplet_to_col(ncol, ncol, ntrpl, Ti, Tj, Tx, Ap, TTi, TTx);
76 :     nz = Ap[ncol]; /* non-zeros in Z'Z crosstab */
77 :     SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nz));
78 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nz));
79 :     Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), TTi, nz);
80 :     Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), TTx, nz);
81 :     Free(Ti); Free(Tj); Free(Tx);
82 :     Free(TTi); Free(TTx);
83 :    
84 :     UNPROTECT(1);
85 :     return val;
86 :     }
87 : bates 73
88 :     extern void ssclme_fill_LIp(int n, const int Parent[], int LIp[]);
89 :    
90 :     SEXP sscCrosstab_L_LI_sizes(SEXP ctab, SEXP permexp)
91 :     {
92 :     SEXP ans = PROTECT(allocVector(INTSXP, 4));
93 :     int *Ai = INTEGER(GET_SLOT(ctab, Matrix_iSym)),
94 :     *Ap = INTEGER(GET_SLOT(ctab, Matrix_pSym)),
95 :     *aa = INTEGER(ans),
96 :     *perm = INTEGER(permexp),
97 :     n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
98 :     *Lp = Calloc(n + 1, int),
99 :     *Parent = Calloc(n, int),
100 :     *Lnz = Calloc(n, int),
101 :     *Flag = Calloc(n, int);
102 :    
103 :     ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag,
104 :     (int *) NULL, (int *) NULL); /* P & Pinv */
105 :     aa[0] = Lp[n];
106 :     ssclme_fill_LIp(n, Parent, Lp);
107 :     aa[1] = Lp[n];
108 :     ssc_symbolic_permute(n, 1, perm, Ap, Ai);
109 :     ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag,
110 :     (int *) NULL, (int *) NULL); /* P & Pinv */
111 :     aa[2] = Lp[n];
112 :     ssclme_fill_LIp(n, Parent, Lp);
113 :     aa[3] = Lp[n];
114 :     Free(Flag); Free(Lnz); Free(Parent); Free(Lp);
115 :     UNPROTECT(1);
116 :     return ans;
117 :     }
118 : bates 85
119 :     /**
120 :     * Determine the inverse of the permutation of the rows of the sparse,
121 :     * column-oriented matrix represented by m, n, Ap and Ai that will
122 :     * concentrate non-zeros in the lower rows.
123 :     *
124 :     * @param m number of rows in the matrix
125 :     * @param n number of columns in the matrix
126 :     * @param Ap column pointers in Ai (modified)
127 :     * @param Ai row indices (modified)
128 : bates 86 * @param perm on return contains the permutation of the rows
129 :     * @param useL when more than one row matches maxrc, use last match
130 : bates 85 */
131 :     static
132 : bates 86 void pair_perm(int m, int n, int Ap[], int Ai[], int iperm[], int useL)
133 : bates 85 {
134 :     int *cc = Calloc(n, int), /* column counts */
135 :     ii, j,
136 :     *rc = Calloc(m, int); /* row counts */
137 : bates 74
138 : bates 85 for (j = 0; j < n; j++) { /* initialize col counts */
139 :     cc[j] = Ap[j+1] - Ap[j];
140 :     }
141 :     for (ii = m - 1; 0 <= ii; ii--) { /* fill iperm from RHS */
142 :     int maxrc, rr;
143 :     int i, mincc, p1, p3;
144 :    
145 :     mincc = m + 1; /* find minimum positive cc */
146 :     for (j = 0; j < n; j++) {
147 :     if (0 < cc[j] && cc[j] < mincc) mincc = cc[j];
148 :     if (mincc < 2) break;
149 :     }
150 :    
151 :     for (i = 0; i < m; i++) rc[i] = 0;
152 :     for (j = 0; j < n; j++) { /* row counts for cols where cc = mincc */
153 :     if (cc[j] == mincc) {
154 :     int p2 = Ap[j+1];
155 :     for (i = Ap[j]; i < p2; i++) rc[Ai[i]]++;
156 :     }
157 :     }
158 :    
159 :     maxrc = -1; /* find last row with rc[i] == max(rc) */
160 :     for (i = 0; i < m; i++) {
161 :     int ic = rc[i];
162 : bates 86 if (ic > maxrc || (useL && ic == maxrc)) {
163 : bates 85 maxrc = ic;
164 :     rr = i;
165 :     }
166 :     }
167 :    
168 : bates 86 iperm[rr] = ii;
169 : bates 85
170 :     p1 = p3 = 0; /* update cc, Ap and Ai */
171 :     for (j = 0; j < n; j++) {
172 :     int p2 = Ap[j+1];
173 :     for (i = Ap[j]; i < p2; i++) {
174 :     if (Ai[i] == rr) {
175 :     cc[j]--;
176 :     } else {
177 :     if (i != p1) Ai[p1] = Ai[i];
178 :     p1++;
179 :     }
180 :     }
181 :     Ap[j] = p3;
182 :     p3 = p1; /* save current pos for next iteration */
183 :     }
184 :     Ap[n] = p3;
185 :     }
186 :     Free(cc); Free(rc);
187 :     }
188 :    
189 : bates 86 SEXP sscCrosstab_groupedPerm(SEXP ctab, SEXP useLast)
190 : bates 74 {
191 :     SEXP
192 :     GpSlot = GET_SLOT(ctab, Matrix_GpSym),
193 :     iSlot = GET_SLOT(ctab, Matrix_iSym),
194 :     pSlot = GET_SLOT(ctab, Matrix_pSym);
195 :     int *Ai = INTEGER(iSlot),
196 :     *Ap = INTEGER(pSlot),
197 :     *Gp = INTEGER(GpSlot),
198 : bates 86 useL = asLogical(useLast),
199 : bates 85 i,
200 : bates 74 n = length(pSlot) - 1, /* number of columns */
201 :     nf = length(GpSlot) - 1, /* number of factors */
202 : bates 85 *np = Calloc(n + 1, int), /* column pointers */
203 :     *ni = Calloc(length(iSlot) - n, int); /* row indices */
204 : bates 74 SEXP ans = PROTECT(allocVector(INTSXP, n));
205 :    
206 :     if (toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L')
207 :     error("Lower triangle required in sscCrosstab object");
208 :    
209 : bates 85 for (i = 0; i < n; i++) {
210 : bates 86 INTEGER(ans)[i] = i; /* initialize permutation to identity */
211 : bates 85 }
212 :     np[0] = 0;
213 : bates 77
214 : bates 85 for (i = 1; i < nf; i++) { /* adjacent pairs of grouping factors */
215 :     int j, k, p0 = 0, p1 = Gp[i-1], p2 = Gp[i], p3 = Gp[i+1];
216 :    
217 :     for (j = p1; j < p2; j++) { /* for this set of columns */
218 :     int lk = Ap[j+1];
219 :     for (k = Ap[j]; k < lk; k++) {
220 :     int ii = Ai[k];
221 :     if (p2 <= ii && ii < p3) { /* check the row */
222 :     ni[p0++] = ii - p2;
223 : bates 77 }
224 : bates 74 }
225 : bates 85 np[j + 1 - p1] = p0;
226 : bates 74 }
227 : bates 86 pair_perm(p3 - p2, p2 - p1, np, ni, INTEGER(ans) + p2, useL);
228 :     for (j = p2; j < p3; j++) INTEGER(ans)[j] += p2;
229 : bates 74 }
230 : bates 86
231 :     Free(np); Free(ni);
232 : bates 74 UNPROTECT(1);
233 :     return ans;
234 :     }

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