SCM

SCM Repository

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

Diff of /pkg/src/sscCrosstab.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 147, Mon May 3 05:26:01 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.147  
changed lines
  Added in v.148

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