SCM Repository

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

Diff of /pkg/src/sscCrosstab.c

revision 146, Sat May 1 20:27:04 2004 UTC revision 148, Mon May 3 12:49:17 2004 UTC
# Line 72  Line 72
72      return val;      return val;
73  }  }
74
extern void ssclme_fill_LIp(int n, const int Parent[], int LIp[]);

75  SEXP sscCrosstab_L_LI_sizes(SEXP ctab, SEXP permexp)  SEXP sscCrosstab_L_LI_sizes(SEXP ctab, SEXP permexp)
76  {  {
77      SEXP ans = PROTECT(allocVector(INTSXP, 4));      SEXP ans = PROTECT(allocVector(INTSXP, 4));
# Line 103  Line 101
101      return ans;      return ans;
102  }  }
103
/**
* Check for a nested series of grouping factors in the sparse,
* symmetric representation of the pairwise cross-tabulations.
*
* @param n size of pairwise cross-tabulation matrix
* @param nf number of groups of columns in pairwise cross-tabulation
* @param upper non-zero if the upper triangle is stored
* @param Ap array of pointers to columns
* @param Ai row indices
* @param Gp array of pointers to groups
*
* @return 0 for non-nested groups, 1 for nested groups
*/
static
int ctab_isNested(int n, int nf, int upper,
const int Ap[], const int Ai[], const int Gp[])
{
if (nf > 1) {  /* single factor always nested */
int  i;
if (upper) {
int *nnz = (int *) R_alloc(n, sizeof(int)), nz = Ap[n];
/* count number of nonzeros in each row */
for (i = 0; i < n; i++) nnz[i] = 0;
for (i = 0; i < nz; i++) nnz[Ai[i]]++;
for (i = 0; i < nf; i++) {
int j, p2 = Gp[i+1], target = nf - i;
for (j = Gp[i]; j < p2; j++) {
if (nnz[j] != target) return 0;
}
}
} else {                /* lower triangle - the easy case */
for (i = 0; i < nf; i++) {
int j, p2 = Gp[i+1], target = nf - i;
for (j = Gp[i]; j < p2; j++) {
if ((Ap[j+1] - Ap[j]) != target)
return 0;
}
}
}
}
return 1;
}

/**
* Determine the inverse of the permutation of the rows of the sparse,
* column-oriented matrix represented by m, n, Ap and Ai that will
* concentrate non-zeros in the lower rows.
*
* @param m number of rows in the matrix
* @param n number of columns in the matrix
* @param Ap column pointers in Ai (modified)
* @param Ai row indices (modified)
* @param rperm on return contains the permutation of the rows
* @param cperm if non-null, on return contains the permutation of the columns
*/
104  static  static
105  void pair_perm(int m, int n, int Ap[], int Ai[], int rperm[], int cperm[])  void col_metis_order(int j0, int j1, int i2,
106                         const int Tp[], const int Ti[], int ans[])
107  {  {
108      int *cc = Calloc(n, int),   /* column counts */      int j, nz = 0;              /* count off-diagonal pairs */
109          *cm = Calloc(n, int),   /* column removed */      for (j = j0; j < j1; j++) { /* columns of interest */
110          *sm = Calloc(m, int),   /* sum of column removals for this row */          int ii, nr = 0, p2 = Tp[j + 1];
111          ii, j, pc,          for (ii = Tp[j]; ii < p2; ii++) {
112          *rc = Calloc(m, int);   /* row counts */              int i = Ti[ii];
113                if (j1 <= i && i < i2) nr++; /* verify row index */
114      for (j = 0; j < n; j++) {   /* initialize col counts */          }
115          cc[j] = Ap[j+1] - Ap[j];          nz += (nr * (nr - 1))/2; /* add number of pairs of rows */
116          cm[j] = 0;      }
117      }      if (nz > 0) {               /* Form an ssc Matrix */
118      pc = 0;          int j, n = i2 - j1,     /* number of rows */
119      for (ii = m - 1; 0 <= ii; ii--) { /* fill rperm from RHS */              nnz = n + nz, pos;
120          int maxrc, rr;          int *Ap = Calloc(n + 1, int),
121          int i, mincc, p1, p3;              *Ai = Calloc(nnz, int),
122                *Tj = Calloc(nnz, int),
123          mincc = m + 1;          /* find minimum positive cc */              *TTi = Calloc(nnz, int);
124          for (j = 0; j < n; j++) {          double                  /* needed for triplet_to_col */
125              if (0 < cc[j] && cc[j] < mincc) mincc = cc[j];              *Ax = Calloc(nnz, double), /* FIXME: change triplet_to_col */
126              if (mincc < 2) break;              *Tx = Calloc(nnz, double); /* to check for null pointers */
127          }          idxtype *perm = Calloc(n, idxtype),
128                *iperm = Calloc(n, idxtype);
129          for (i = 0; i < m; i++) {sm[i] = rc[i] = 0;}
130          for (j = 0; j < n; j++) { /* row counts for cols where cc = mincc */          for (j = 0; j < n; j++) { /* diagonals */
131              if (cc[j] == mincc) {              TTi[j] = Tj[j] = j;
132                  int p2 = Ap[j+1];              Tx[j] = 1.;
133                  for (i = Ap[j]; i < p2; i++) {          }
134                      rc[Ai[i]]++;          pos = n;
135                      sm[Ai[i]] += cm[j];          for (j = j0; j < j1; j++) { /* create the pairs */
136                  }              int ii, nr = 0, p2 = Tp[j + 1];
137                for (ii = Tp[j]; ii < p2; ii++) {
138                    int r1 = Ti[ii], i1;
139                    if (j1 <= r1 && r1 < i2) {
140                        for (i1 = ii + 1; i1 < p2; i1++) {
141                            int r2 = Ti[i1];
142                            if (r2 < i2) {
143                                TTi[pos] = r2 - j1;
144                                Tj[pos] = r1 - j1;
145                                Tx[pos] = 1.;
146                                pos++;
147              }              }
148          }          }
maxrc = -1;             /* find rows with rc[i] == max(rc) */
for (i = 0; i < m; i++) {
int ic = rc[i];
/* Choose first on row count.  Ties go
* to smaller sum of moved.  Ties
* there go to the last one. */
if (ic > maxrc || (ic == maxrc && sm[i] >= sm[rr])) {
maxrc = ic;
rr = i;
}
}

rperm[rr] = ii;

p1 = p3 = 0;            /* update cc, Ap and Ai */
for (j = 0; j < n; j++) {
int p2 = Ap[j+1];
for (i = Ap[j]; i < p2; i++) {
if (Ai[i] == rr) {
cc[j]--; cm[j]++; /* move from count to removed */
if (cperm && cc[j] < 1) cperm[j] = pc++;
} else {
if (i != p1) Ai[p1] = Ai[i];
p1++;
149                  }                  }
150              }              }
Ap[j] = p3;
p3 = p1;            /* save current pos for next iteration */
151          }          }
152          Ap[n] = p3;          triplet_to_col(n, n, nnz, TTi, Tj, Tx, Ap, Ai, Ax);
153            ssc_metis_order(n, Ap, Ai, perm, iperm);
154            for (j = j1; j < i2; j++) ans[j] = j1 + iperm[j - j1];
155            Free(Tx); Free(Ax); Free(TTi); Free(Tj); Free(Ai); Free(Ap);
156            Free(perm); Free(iperm);
157      }      }
Free(cc); Free(cm); Free(rc);
158  }  }
159
160  SEXP sscCrosstab_groupedPerm(SEXP ctab)  SEXP sscCrosstab_groupedPerm(SEXP ctab)

Legend:
 Removed from v.146 changed lines Added in v.148