# SCM Repository

[matrix] Diff of /pkg/src/sscCrosstab.c
 [matrix] / pkg / src / sscCrosstab.c

# Diff of /pkg/src/sscCrosstab.c

revision 83, Thu Apr 15 02:51:07 2004 UTC revision 85, Thu Apr 15 20:50:57 2004 UTC
# Line 116  Line 116
116      return ans;      return ans;
117  }  }
118
119  /* FIXME:  /**
120   *  Split main part into a function that takes nr, nc, p, i, and    * Determine the inverse of the permutation of the rows of the sparse,
121   *  returns perm of size nr.    * column-oriented matrix represented by m, n, Ap and Ai that will
122   *  Call function for successive pairs.    * concentrate non-zeros in the lower rows.
123   *  Invert the permutation before return.    *
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
137        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  SEXP sscCrosstab_groupedPerm(SEXP ctab)  SEXP sscCrosstab_groupedPerm(SEXP ctab)
189  {  {
190      SEXP      SEXP
# Line 131  Line 194
194      int *Ai = INTEGER(iSlot),      int *Ai = INTEGER(iSlot),
195          *Ap = INTEGER(pSlot),          *Ap = INTEGER(pSlot),
196          *Gp = INTEGER(GpSlot),          *Gp = INTEGER(GpSlot),
197          nl1 = Gp[1],            /* number of levels of first factor */          i,
*icounts,               /* counts for rows */
*jcounts = Calloc(nl1, int), /* counts for columns */
nl2,                    /* delay initializing until nf is known */
ii, j, jj,
198          n = length(pSlot) - 1,  /* number of columns */          n = length(pSlot) - 1,  /* number of columns */
199          nf = length(GpSlot) - 1, /* number of factors */          nf = length(GpSlot) - 1, /* number of factors */
200          nz = length(iSlot),     /* number of non-zeros */          *np = Calloc(n + 1, int), /* column pointers */
201          nod = nz - n,           /* number of off-diagonals */          *ni = Calloc(length(iSlot) - n, int); /* row indices */
*np = Calloc(nl1 + 1, int), /* new array pointers */
*ni = Calloc(nod, int), /* new array row indices */
p1, p3;
202      SEXP ans = PROTECT(allocVector(INTSXP, n));      SEXP ans = PROTECT(allocVector(INTSXP, n));
203      int *perm = INTEGER(ans);      int *iperm = Calloc(n, int);
204
205      if (toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L')      if (toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L')
206          error("Lower triangle required in sscCrosstab object");          error("Lower triangle required in sscCrosstab object");
207
208      for (j = 0; j < n; j++) perm[j] = j; /* initialize permutation to identity */      for (i = 0; i < n; i++) {
209            iperm[i] = i;    /* initialize inverse permutation to identity */
if (nf > 1) {
nl2 = Gp[2] - nl1;      /* number of levels of factor 2 */
icounts = Calloc(nl2, int);
np[0] = 0;
for (j = 0; j < nl1; j++) jcounts[j] = 0;
p1 = 0;                 /* copy off-diagonals from first nl1 cols */
for (j = 0; j < nl1; j++) { /* and 0 <= row - nl1 < nl2 */
int p3 = Ap[j + 1];
for (jj = Ap[j]; jj < p3; jj++) {
int i = Ai[jj] - nl1;
if (0 <= i && i < nl2) {
ni[p1++] = i;
jcounts[j]++;
}
}
np[j+1] = p1;
210          }          }
211        np[0] = 0;
212
213          for (ii = 1; ii <= nl2; ii++) {      for (i = 1; i < nf; i++) {  /* adjacent pairs of grouping factors */
214              int maxrc, rr;          int j, k, p0 = 0, p1 = Gp[i-1], p2 = Gp[i], p3 = Gp[i+1];
215              int i, minjc = nl2+1;       /* find minimum positive jcount */
216              for (j = 0; j < nl1; j++) {          for (j = p1; j < p2; j++) { /* for this set of columns */
217                  if (jcounts[j] > 0 && jcounts[j] < minjc) minjc = jcounts[j];              int lk = Ap[j+1];
218                  if (minjc == 1) break;              for (k = Ap[j]; k < lk; k++) {
219              }                  int ii = Ai[k];
220                                  /* accumulate the row counts on cols where jcount=minjc */                  if (p2 <= ii && ii < p3) { /* check the row */
221              for (i = 0; i < nl2; i++) icounts[i] = 0;                      ni[p0++] = ii - p2;
for (j = 0; j < nl1; j++) {
if (jcounts[j] == minjc) {
int p2 = np[j+1];
for (i = np[j]; i < p2; i++) icounts[ni[i]]++;
}
}
/* find the last row with count==max(icount) */
maxrc = -1;
for (i = 0; i < nl2; i++) {
int ic = icounts[i];
if (ic >= maxrc) {
maxrc = ic;
rr = i;
}
}
perm[nl1 + nl2 - ii] = rr + nl1;
/* update icounts, np and ni */
p1 = p3 = 0;
for (j = 0; j < nl1; j++) {
int p2 = np[j+1];
for (i = np[j]; i < p2; i++) {
if (ni[i] == rr) {
jcounts[j]--;
} else {
if (i != p1) ni[p1] = ni[i];
p1++;
222                      }                      }
223                  }                  }
224                  np[j] = p3;              np[j + 1 - p1] = p0;
p3 = p1;        /* save the count for the next iteration */
225              }              }
226              np[nl1] = p3;          pair_perm(p3 - p2, p2 - p1, np, ni, iperm + p2);
227            for (j = p2; j < p3; j++) iperm[j] += p2;
228          }          }
229          Free(icounts);
230        for(i = 0; i < n; i++) {    /* invert the permutation */
231            INTEGER(ans)[iperm[i]] = i;
232      }      }
233      Free(np); Free(ni); Free(jcounts);      Free(np); Free(ni); Free(iperm);
234      UNPROTECT(1);      UNPROTECT(1);
235      return ans;      return ans;
236  }  }

Legend:
 Removed from v.83 changed lines Added in v.85