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 85 - (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 :     * @param iperm on return contains the inverse permutation
129 :     */
130 :     static
131 :     void pair_perm(int m, int n, int Ap[], int Ai[], int iperm[])
132 :     {
133 :     int *cc = Calloc(n, int), /* column counts */
134 :     ii, j,
135 :     *rc = Calloc(m, int); /* row counts */
136 : bates 74
137 : bates 85 for (j = 0; j < n; j++) { /* initialize col counts */
138 :     cc[j] = Ap[j+1] - Ap[j];
139 :     }
140 :     for (ii = m - 1; 0 <= ii; ii--) { /* fill iperm from RHS */
141 :     int maxrc, rr;
142 :     int i, mincc, p1, p3;
143 :    
144 :     mincc = m + 1; /* find minimum positive cc */
145 :     for (j = 0; j < n; j++) {
146 :     if (0 < cc[j] && cc[j] < mincc) mincc = cc[j];
147 :     if (mincc < 2) break;
148 :     }
149 :    
150 :     for (i = 0; i < m; i++) rc[i] = 0;
151 :     for (j = 0; j < n; j++) { /* row counts for cols where cc = mincc */
152 :     if (cc[j] == mincc) {
153 :     int p2 = Ap[j+1];
154 :     for (i = Ap[j]; i < p2; i++) rc[Ai[i]]++;
155 :     }
156 :     }
157 :    
158 :     maxrc = -1; /* find last row with rc[i] == max(rc) */
159 :     for (i = 0; i < m; i++) {
160 :     int ic = rc[i];
161 :     if (ic >= maxrc) {
162 :     maxrc = ic;
163 :     rr = i;
164 :     }
165 :     }
166 :    
167 :     iperm[ii] = rr;
168 :    
169 :     p1 = p3 = 0; /* update cc, Ap and Ai */
170 :     for (j = 0; j < n; j++) {
171 :     int p2 = Ap[j+1];
172 :     for (i = Ap[j]; i < p2; i++) {
173 :     if (Ai[i] == rr) {
174 :     cc[j]--;
175 :     } else {
176 :     if (i != p1) Ai[p1] = Ai[i];
177 :     p1++;
178 :     }
179 :     }
180 :     Ap[j] = p3;
181 :     p3 = p1; /* save current pos for next iteration */
182 :     }
183 :     Ap[n] = p3;
184 :     }
185 :     Free(cc); Free(rc);
186 :     }
187 :    
188 : bates 74 SEXP sscCrosstab_groupedPerm(SEXP ctab)
189 :     {
190 :     SEXP
191 :     GpSlot = GET_SLOT(ctab, Matrix_GpSym),
192 :     iSlot = GET_SLOT(ctab, Matrix_iSym),
193 :     pSlot = GET_SLOT(ctab, Matrix_pSym);
194 :     int *Ai = INTEGER(iSlot),
195 :     *Ap = INTEGER(pSlot),
196 :     *Gp = INTEGER(GpSlot),
197 : bates 85 i,
198 : bates 74 n = length(pSlot) - 1, /* number of columns */
199 :     nf = length(GpSlot) - 1, /* number of factors */
200 : bates 85 *np = Calloc(n + 1, int), /* column pointers */
201 :     *ni = Calloc(length(iSlot) - n, int); /* row indices */
202 : bates 74 SEXP ans = PROTECT(allocVector(INTSXP, n));
203 : bates 85 int *iperm = Calloc(n, int);
204 : bates 74
205 :     if (toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L')
206 :     error("Lower triangle required in sscCrosstab object");
207 :    
208 : bates 85 for (i = 0; i < n; i++) {
209 :     iperm[i] = i; /* initialize inverse permutation to identity */
210 :     }
211 :     np[0] = 0;
212 : bates 77
213 : bates 85 for (i = 1; i < nf; i++) { /* adjacent pairs of grouping factors */
214 :     int j, k, p0 = 0, p1 = Gp[i-1], p2 = Gp[i], p3 = Gp[i+1];
215 :    
216 :     for (j = p1; j < p2; j++) { /* for this set of columns */
217 :     int lk = Ap[j+1];
218 :     for (k = Ap[j]; k < lk; k++) {
219 :     int ii = Ai[k];
220 :     if (p2 <= ii && ii < p3) { /* check the row */
221 :     ni[p0++] = ii - p2;
222 : bates 77 }
223 : bates 74 }
224 : bates 85 np[j + 1 - p1] = p0;
225 : bates 74 }
226 : bates 85 pair_perm(p3 - p2, p2 - p1, np, ni, iperm + p2);
227 :     for (j = p2; j < p3; j++) iperm[j] += p2;
228 : bates 74 }
229 : bates 85
230 :     for(i = 0; i < n; i++) { /* invert the permutation */
231 :     INTEGER(ans)[iperm[i]] = i;
232 :     }
233 :     Free(np); Free(ni); Free(iperm);
234 : bates 74 UNPROTECT(1);
235 :     return ans;
236 :     }

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