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 73, Tue Apr 13 15:53:38 2004 UTC revision 74, Wed Apr 14 16:37:24 2004 UTC
# Line 44  Line 44 
44      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
45      dims = INTEGER(GET_SLOT(val, Matrix_DimSym));      dims = INTEGER(GET_SLOT(val, Matrix_DimSym));
46      dims[0] = dims[1] = ncol;      dims[0] = dims[1] = ncol;
     SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar("U")));  
47      ntrpl = nfc2 * nobs + ncol;      ntrpl = nfc2 * nobs + ncol;
48      Ti = Calloc(ntrpl, int); Tj = Calloc(ntrpl, int); TTi = Calloc(ntrpl, int);      Ti = Calloc(ntrpl, int); Tj = Calloc(ntrpl, int); TTi = Calloc(ntrpl, int);
49      Tx = Calloc(ntrpl, double); TTx = Calloc(ntrpl, double);      Tx = Calloc(ntrpl, double); TTx = Calloc(ntrpl, double);
# Line 116  Line 115 
115      UNPROTECT(1);      UNPROTECT(1);
116      return ans;      return ans;
117  }  }
118    
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    }

Legend:
Removed from v.73  
changed lines
  Added in v.74

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