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 82 - (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 74
119 :     SEXP sscCrosstab_groupedPerm(SEXP ctab)
120 :     {
121 :     SEXP
122 :     GpSlot = GET_SLOT(ctab, Matrix_GpSym),
123 :     iSlot = GET_SLOT(ctab, Matrix_iSym),
124 :     pSlot = GET_SLOT(ctab, Matrix_pSym);
125 :     int *Ai = INTEGER(iSlot),
126 :     *Ap = INTEGER(pSlot),
127 :     *Gp = INTEGER(GpSlot),
128 :     nl1 = Gp[1], /* number of levels of first factor */
129 : bates 77 *icounts, /* counts for rows */
130 :     *jcounts = Calloc(nl1, int), /* counts for columns */
131 :     nl2, /* delay initializing until nf is known */
132 :     ii, j, jj,
133 : bates 74 n = length(pSlot) - 1, /* number of columns */
134 :     nf = length(GpSlot) - 1, /* number of factors */
135 :     nz = length(iSlot), /* number of non-zeros */
136 :     nod = nz - n, /* number of off-diagonals */
137 : bates 77 *np = Calloc(nl1 + 1, int), /* new array pointers */
138 :     *ni = Calloc(nod, int), /* new array row indices */
139 :     p1, p3;
140 : bates 74 SEXP ans = PROTECT(allocVector(INTSXP, n));
141 :     int *perm = INTEGER(ans);
142 :    
143 :     if (toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L')
144 :     error("Lower triangle required in sscCrosstab object");
145 :    
146 : bates 77 for (j = 0; j < n; j++) perm[j] = j; /* initialize permutation to identity */
147 :    
148 :     if (nf > 1) {
149 : bates 82 nl2 = Gp[2] - nl1; /* number of levels of factor 2 */
150 :     icounts = Calloc(nl2, int);
151 : bates 77 np[0] = 0;
152 :     for (j = 0; j < nl1; j++) jcounts[j] = 0;
153 :     p1 = 0; /* copy off-diagonals from first nl1 cols */
154 : bates 82 for (j = 0; j < nl1; j++) { /* and 0 <= row - nl1 < nl2 */
155 : bates 77 int p3 = Ap[j + 1];
156 :     for (jj = Ap[j]; jj < p3; jj++) {
157 : bates 82 int i = Ai[jj] - nl1;
158 :     if (0 <= i && i < nl2) {
159 :     ni[p1++] = i;
160 : bates 77 jcounts[j]++;
161 :     }
162 : bates 74 }
163 : bates 77 np[j+1] = p1;
164 : bates 74 }
165 :    
166 : bates 77 for (ii = 1; ii <= nl2; ii++) {
167 :     int maxrc, rr;
168 :     int i, minjc = nl2+1; /* find minimum positive jcount */
169 :     for (j = 0; j < nl1; j++) {
170 :     if (jcounts[i] > 0 && jcounts[j] < minjc) minjc = jcounts[j];
171 :     }
172 :     /* accumulate the row counts on cols where jcount=minjc */
173 :     for (i = 0; i < nl2; i++) icounts[i] = 0;
174 :     for (j = 0; j < nl1; j++) {
175 :     if (jcounts[j] == minjc) {
176 :     int p2 = np[j+1];
177 :     for (i = np[j]; i < p2; i++) icounts[ni[i]]++;
178 : bates 74 }
179 :     }
180 : bates 82 /* find the last row with count==max(icount) */
181 : bates 77 maxrc = -1;
182 :     for (i = 0; i < nl2; i++) {
183 :     int ic = icounts[i];
184 :     if (ic >= maxrc) {
185 :     maxrc = ic;
186 :     rr = i;
187 :     }
188 : bates 74 }
189 : bates 82 perm[nl1 + nl2 - ii] = rr + nl1;
190 : bates 77 /* update icounts, np and ni */
191 :     p1 = p3 = 0;
192 :     for (j = 0; j < nl1; j++) {
193 :     int p2 = np[j+1];
194 :     for (i = np[j]; i < p2; i++) {
195 :     if (ni[i] != rr) {
196 :     if (i != p1) ni[p1] = ni[i];
197 :     p1++;
198 :     }
199 :     }
200 :     np[j] = p3;
201 : bates 82 p3 = p1; /* save the count for the next iteration */
202 : bates 74 }
203 :     }
204 : bates 77 Free(icounts);
205 : bates 74 }
206 : bates 77 Free(np); Free(ni); Free(jcounts);
207 : bates 74 UNPROTECT(1);
208 :     return ans;
209 :     }

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