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 74 - (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 :     static
120 :     void make_icounts(int nod, int ni, const int i[], const int j[],
121 :     const char ind[], int icounts[])
122 :     {
123 :     int k, *lastj = memset(Calloc(ni, int), 0, sizeof(int) * ni);
124 :    
125 :     memset(icounts, 0, sizeof(int) * ni);
126 :     for (k = 0; k < nod; k++) {
127 :     int ik = i[k], jk = j[k];
128 :     if (ind[k] && jk != lastj[ik]) {
129 :     icounts[ik]++;
130 :     lastj[ik] = jk;
131 :     }
132 :     }
133 :     Free(lastj);
134 :     }
135 :    
136 :     SEXP sscCrosstab_groupedPerm(SEXP ctab)
137 :     {
138 :     SEXP
139 :     GpSlot = GET_SLOT(ctab, Matrix_GpSym),
140 :     iSlot = GET_SLOT(ctab, Matrix_iSym),
141 :     pSlot = GET_SLOT(ctab, Matrix_pSym);
142 :     int *Ai = INTEGER(iSlot),
143 :     *Ap = INTEGER(pSlot),
144 :     *Gp = INTEGER(GpSlot),
145 :     nl1 = Gp[1], /* number of levels of first factor */
146 :     *icounts = Calloc(nl1, int),
147 :     nl2, *jcounts,
148 :     j, jj,
149 :     n = length(pSlot) - 1, /* number of columns */
150 :     nf = length(GpSlot) - 1, /* number of factors */
151 :     nz = length(iSlot), /* number of non-zeros */
152 :     nod = nz - n, /* number of off-diagonals */
153 :     nuse = nod,
154 :     *iv = Calloc(nod, int),
155 :     *jv = Calloc(nod, int),
156 :     p1, p2;
157 :     char *active = (char *) memset(Calloc(n, char), 1, (size_t) n),
158 :     *ind = (char *) memset(Calloc(nod, char), 1, (size_t) nod);
159 :    
160 :     SEXP ans = PROTECT(allocVector(INTSXP, n));
161 :     int *perm = INTEGER(ans);
162 :    
163 :     if (toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L')
164 :     error("Lower triangle required in sscCrosstab object");
165 :     if (nf < 2) error("At least two factors are required");
166 :     nl2 = Gp[2] - Gp[1];
167 :     jcounts = Calloc(nl2, int);
168 :    
169 :     p1 = 0; p2 = nl1; /* copy and expand off-diagonal indices */
170 :     for (j = 0; j < p2; j++) {
171 :     int p3 = Ap[j + 1];
172 :     for (jj = Ap[j]; jj < p3; jj++) {
173 :     int i = Ai[jj];
174 :     /* FIXME: iv and jv are named backwards. Reconcile this. */
175 :     if (i > j) {
176 :     jv[p1] = i;
177 :     iv[p1] = j;
178 :     p1++;
179 :     }
180 :     }
181 :     }
182 :     if (p1 != nod) error("Logic error - counts of off-diagonals disagree");
183 :    
184 :     p1 = 0; p2 = Gp[2] - 1; /* fill 1st level from LHS, 2nd from RHS */
185 :     while (nuse > 0) {
186 :     int k, mcount, mind;
187 :     make_icounts(nod, nl1, iv, jv, ind, icounts);
188 :    
189 :     mind = -1;
190 :     mcount = Gp[2]; /* must be bigger than Gp[2] - Gp[1] */
191 :     for (k = 0; k < nl1; k++) { /* determine column with min count */
192 :     int ic = icounts[k];
193 :     if (active[k]) {
194 :     if (ic < 1) { /* remove active columns with zero counts */
195 :     perm[p1++] = k;
196 :     active[k] = 0;
197 :     } else if (ic < mcount) {
198 :     mcount = ic;
199 :     mind = k;
200 :     }
201 :     }
202 :     }
203 :     if (mind < 0)
204 :     error("Logic error - ran out of columns before nuse == 0");
205 :    
206 :     /* Count rows for j */
207 :     memset(jcounts, 0, sizeof(int) * nl2);
208 :     for (k = 0; k < nod; k++) {
209 :     if (ind[k] && icounts[iv[k]] == mcount) {
210 :     jcounts[jv[k] - nl1]++;
211 :     }
212 :     }
213 :    
214 :     mcount = -1; mind = -1; /* determine j with max count */
215 :     for (k = 0; k < nl2; k++) {
216 :     int jc = jcounts[k];
217 :     /* use >= below so ties give last row */
218 :     if (active[k + nl1] && jc >= mcount) {
219 :     mcount = jc;
220 :     mind = k;
221 :     }
222 :     }
223 :     if (mind < 0)
224 :     error("Logic error - ran out of rows before nuse == 0");
225 :    
226 :     perm[p2--] = mind + nl1;
227 :     for (k = 0; k < nod; k++) {
228 :     if (jv[k] - nl1 == mind) {
229 :     ind[k] = 0;
230 :     nuse--;
231 :     }
232 :     }
233 :     }
234 :    
235 :     Free(ind); Free(active); Free(jv); Free(iv); Free(jcounts); Free(icounts);
236 :     UNPROTECT(1);
237 :     return ans;
238 :     }

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