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 127, Sun Apr 25 18:06:27 2004 UTC revision 143, Fri Apr 30 20:33:01 2004 UTC
# Line 117  Line 117 
117  }  }
118    
119  /**  /**
120     * Check for a nested series of grouping factors in the sparse,
121     * symmetric representation of the pairwise cross-tabulations.
122     *
123     * @param n size of pairwise cross-tabulation matrix
124     * @param nf number of groups of columns in pairwise cross-tabulation
125     * @param upper non-zero if the upper triangle is stored
126     * @param Ap array of pointers to columns
127     * @param Ai row indices
128     * @param Gp array of pointers to groups
129     *
130     * @return 0 for non-nested groups, 1 for nested groups
131     */
132    static
133    int ctab_isNested(int n, int nf, int upper,
134                      const int Ap[], const int Ai[], const int Gp[])
135    {
136        if (nf > 1) {  /* single factor always nested */
137            int  i;
138            if (upper) {
139                int *nnz = (int *) R_alloc(n, sizeof(int)), nz = Ap[n];
140                                    /* count number of nonzeros in each row */
141                for (i = 0; i < n; i++) nnz[i] = 0;
142                for (i = 0; i < nz; i++) nnz[Ai[i]]++;
143                for (i = 0; i < nf; i++) {
144                    int j, p2 = Gp[i+1], target = nf - i;
145                    for (j = Gp[i]; j < p2; j++) {
146                        if (nnz[j] != target) return 0;
147                    }
148                }
149            } else {                /* lower triangle - the easy case */
150                for (i = 0; i < nf; i++) {
151                    int j, p2 = Gp[i+1], target = nf - i;
152                    for (j = Gp[i]; j < p2; j++) {
153                        if ((Ap[j+1] - Ap[j]) != target)
154                            return 0;
155                    }
156                }
157            }
158        }
159        return 1;
160    }
161    
162    /**
163    * Determine the inverse of the permutation of the rows of the sparse,    * Determine the inverse of the permutation of the rows of the sparse,
164    * column-oriented matrix represented by m, n, Ap and Ai that will    * column-oriented matrix represented by m, n, Ap and Ai that will
165    * concentrate non-zeros in the lower rows.    * concentrate non-zeros in the lower rows.
# Line 208  Line 251 
251          i,          i,
252          n = length(pSlot) - 1,  /* number of columns */          n = length(pSlot) - 1,  /* number of columns */
253          nf = length(GpSlot) - 1, /* number of factors */          nf = length(GpSlot) - 1, /* number of factors */
254          *np = Calloc(n + 1, int), /* column pointers */          up;
         *ni = Calloc(length(iSlot) - n, int); /* row indices */  
255      SEXP ans = PROTECT(allocVector(INTSXP, n));      SEXP ans = PROTECT(allocVector(INTSXP, n));
256    
257      if (toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L')      up = toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L';
258          error("Lower triangle required in sscCrosstab object");      if (nf > 1 && up) {                 /* transpose */
259            int nz = length(iSlot);
260            int *ai = Calloc(nz, int),
261                *ap = Calloc(n + 1, int);
262            double *ax = Calloc(nz, double);
263    
264            csc_components_transpose(n, n, nz, Ap, Ai,
265                                     REAL(GET_SLOT(ctab, Matrix_xSym)),
266                                     ap, ai, ax);
267            Ap = ap;
268            Ai = ai;
269            Free(ax);               /* don't need values, only positions */
270        }
271      for (i = 0; i < n; i++) {      for (i = 0; i < n; i++) {
272          INTEGER(ans)[i] = i;    /* initialize permutation to identity */          INTEGER(ans)[i] = i;    /* initialize permutation to identity */
273      }      }
274      np[0] = 0;      if (!ctab_isNested(n, nf, 0, Ap, Ai, Gp)) {
275            int *np = Calloc(n + 1, int), /* column pointers */
276                *ni = Calloc(length(iSlot) - n, int); /* row indices */
277    
278            np[0] = 0;
279      for (i = 1; i < nf; i++) {  /* adjacent pairs of grouping factors */      for (i = 1; i < nf; i++) {  /* adjacent pairs of grouping factors */
280          int j, k, p0 = 0, p1 = Gp[i-1], p2 = Gp[i], p3 = Gp[i+1];          int j, k, p0 = 0, p1 = Gp[i-1], p2 = Gp[i], p3 = Gp[i+1];
281    
# Line 234  Line 290 
290              np[j + 1 - p1] = p0;              np[j + 1 - p1] = p0;
291          }          }
292          pair_perm(p3 - p2, p2 - p1, np, ni,          pair_perm(p3 - p2, p2 - p1, np, ni,
293                    INTEGER(ans) + p2, INTEGER(ans));                        INTEGER(ans) + p2,
294                          (i == 1) ? INTEGER(ans) + p1 : (int *) 0);
295          for (j = p2; j < p3; j++) INTEGER(ans)[j] += p2;          for (j = p2; j < p3; j++) INTEGER(ans)[j] += p2;
296      }      }
   
297      Free(np); Free(ni);      Free(np); Free(ni);
298        }
299        if (nf > 1 && up) {Free(Ap); Free(Ai);}
300      UNPROTECT(1);      UNPROTECT(1);
301      return ans;      return ans;
302  }  }

Legend:
Removed from v.127  
changed lines
  Added in v.143

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