# SCM Repository

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

# Diff of /pkg/src/bCrosstab.c

revision 347, Tue Nov 23 15:37:58 2004 UTC revision 356, Sat Nov 27 01:45:39 2004 UTC
# Line 1  Line 1
1    /* TODO: 1) combine diag_update and offdiag_update - DONE
2             2) Move Parent_inverse from tscMatrix to here
3             3) Modify the creation of the L diagonals to include a call
4             to ldl_numeric to get Li.
5    */
6
7  #include "bCrosstab.h"  #include "bCrosstab.h"
#include "ldl.h"
8
9  /**  /**
10   * Calculate the zero-based index in a row-wise packed lower   * Calculate the zero-based index in a row-wise packed lower
# Line 25  Line 30
30   * @param perm 0-based permutation vector of length max(i)   * @param perm 0-based permutation vector of length max(i)
31   */   */
32  static R_INLINE  static R_INLINE
33  void ind_permute(int i[], int nnz, int perm[])  void ind_permute(int i[], int nnz, const int perm[])
34  {  {
35      int j;      int j;
36      for (j = 0; j < nnz; j++) i[j] = perm[i[j]];      for (j = 0; j < nnz; j++) i[j] = perm[i[j]];
37  }  }
38
39  /**  /**
40   * Force indices to be in the lower triangle of a matrix   * Force indices to be in the upper triangle of a matrix
41   *   *
42   * @param i vector of 0-based row indices   * @param i vector of 0-based row indices
43   * @param j vector of 0-based column indices   * @param j vector of 0-based column indices
44   * @param nnz length of index vectors   * @param nnz length of index vectors
45   */   */
46  static R_INLINE  static R_INLINE
47  void make_lower_triangular(int i[], int j[], int nnz)  void make_upper_triangular(int i[], int j[], int nnz)
48  {  {
49      int k;      int k;
50      for (k = 0; k < nnz; k++) {      for (k = 0; k < nnz; k++) {
51          if (i[k] < j[k]) {          if (i[k] > j[k]) {
52              int tmp = i[k];              int tmp = i[k];
53              i[k] = j[k];              i[k] = j[k];
54              j[k] = tmp;              j[k] = tmp;
# Line 63  Line 68
68  SEXP make_named_list(int n, char **names)  SEXP make_named_list(int n, char **names)
69  {  {
70      SEXP ans = PROTECT(allocVector(VECSXP, n)),      SEXP ans = PROTECT(allocVector(VECSXP, n)),
71          nms = PROTECT(allocVector(CHARSXP, n));          nms = PROTECT(allocVector(STRSXP, n));
72      int i;      int i;
73
74      for (i = 0; i < n; i++) nms[i] = mkChar(names[i]);      for (i = 0; i < n; i++) SET_STRING_ELT(nms, i, mkChar(names[i]));
75      setAttrib(ans, R_NamesSymbol, nms);      setAttrib(ans, R_NamesSymbol, nms);
76      UNPROTECT(2);      UNPROTECT(2);
77      return ans;      return ans;
78  }  }
79
80    static R_INLINE
81    int check_csc_index(const int p[], const int i[], int row, int col)
82    {
83        int k, k2 = p[col + 1];
84                                    /* use a linear search for now */
85                                    /* perhaps replace by bsearch later */
86        for (k = p[col]; k < k2; k++) {
87            if (i[k] == row) return k;
88        }
89        return -1;
90    }
91
92  /**  /**
93   * Update a diagonal block   * Update a diagonal block
94   *   *
95   * @param ctab pointer to a blocked crosstabulation object   * @param ctab pointer to a blocked crosstabulation object
96   * @param jm1 index of updating column   * @param j index of updating column
97   * @param i index of diagonal block to be updated   * @param i index of diagonal block to be updated
98   */   */
99  static void diag_update(ctab, jm1, i)  static void diag_update(SEXP ctab, int j, int i)
100  {  {
101      SEXP db = VECTOR_ELT(ctab, Lind(i, i)),      SEXP db = VECTOR_ELT(ctab, Lind(i, i)),
102          jb = VECTOR_ELT(ctab, Lind(i, jm1));          jb = VECTOR_ELT(ctab, Lind(i, j));
103      SEXP dpp = GET_SLOT(db, Matrix_pSym),      SEXP dpp = GET_SLOT(db, Matrix_pSym),
104          jpp = GET_SLOT(jb, Matrix_pSym);          jpp = GET_SLOT(jb, Matrix_pSym);
105      int nci = length(dpp) - 1,      int *di = INTEGER(GET_SLOT(db, Matrix_iSym)),
ncj = length(jpp) - 1,
106          *dp = INTEGER(dpp),          *dp = INTEGER(dpp),
107          *di = INTEGER(GET_SLOT(db), Matrix_iSym),          *ji = INTEGER(GET_SLOT(jb, Matrix_iSym)),
108          *jp = INTEGER(jpp),          *jp = INTEGER(jpp),
109          *ji = INTEGER(GET_SLOT(jb), Matrix_iSym);          dnc = length(dpp) - 1,
110            jnc = length(jpp) - 1;
111        int jj, extra;
112                                    /* bound the number of extra elements */
113        extra = 0;
114        for (jj = 0; jj < jnc; jj++) {
115            int k, kk, k2 = jp[jj + 1];
116            for (k = jp[jj]; k < k2; k++) {
117                for (kk = k; kk < k2; kk++) {
118                    if (check_csc_index(dp, di, ji[k], ji[kk]) < 0) extra++;
119                }
120            }
121        }
122        if (!extra) return;
123        {
124            int pos, nnz = dp[dnc];
125            int ntot = nnz + extra;
126            int *Ai = Calloc(ntot, int),
127                *Ti = Calloc(ntot, int),
128                *Tj = Calloc(ntot, int);
129            double *Ax;
130
131            Memcpy(Ti, di, nnz);    /* make a copy of the row indices */
132            pos = 0;                /* fill in the column indices */
133            for (jj = 0; jj < dnc; jj++) {
134                int j2 = dp[jj + 1];
135                while (pos < j2) {
136                    Tj[pos] = jj;
137                    pos++;
138                }
139            }
140                                    /* add the extra elements */
141            for (jj = 0; jj < jnc; jj++) {
142                int k, kk, k2 = jp[jj + 1];
143                for (k = jp[jj]; k < k2; k++) {
144                    for (kk = k; kk < k2; kk++) {
145                        if (check_csc_index(dp, di, ji[k], ji[kk]) < 0) {
146                            Ti[pos] = ji[k];
147                            Tj[pos] = ji[kk];
148                            pos++;
149                        }
150                    }
151                }
152            }
153            triplet_to_col(dnc, dnc, ntot, Ti, Tj, (double *) NULL,
154                           dp, Ai, (double *) NULL);
155            nnz = dp[dnc];
156            SET_SLOT(db, Matrix_iSym, allocVector(INTSXP, nnz));
157            Memcpy(INTEGER(GET_SLOT(db, Matrix_iSym)), Ai, nnz);
158            SET_SLOT(db, Matrix_xSym, allocVector(REALSXP, nnz));
159            Ax = REAL(GET_SLOT(db, Matrix_xSym));
160            for (j = 0; j < nnz; j++) Ax[j] = 1.;
161            Free(Ai); Free(Ti); Free(Tj);
162            return;
163        }
164  }  }
165
166  /**  /**
167   * Project a column of a pairwise crosstabulation onto the remaining columns   * Update a block
168   *   *
169   * @param ctab pointer to a blocked crosstabulation object   * @param ctab pointer to a blocked crosstabulation object
170   * @param i index of column to project   * @param j index of updating column
171   *   * @param k column index of block to be updated
172   * @return a list of the projected ctab, the permutation of factor   * @param i row index of block to be updated (j < k <= i)
* j+1, and the j+1 diagonal element of L-inverse
173   */   */
174  SEXP bCrosstab_project(SEXP ctab, SEXP jp)  static void block_update(SEXP ctab, int j, int k, int i)
175  {  {
176      char anms[][] = {"ctab", "perm", "Linv"};      SEXP tb = VECTOR_ELT(ctab, Lind(i, k)),
177      SEXP ans = PROTECT(make_named_list(3, anms));          ib = VECTOR_ELT(ctab, Lind(i, j)),
178      int ctbl = length(ctab), i, j = asInteger(jp), jm1, k;          kb = VECTOR_ELT(ctab, Lind(k, j));
179      int nf = ((int) (-1 + sqrt((double)(1 + 8*ctbl))))/2;      SEXP tpp = GET_SLOT(tb, Matrix_pSym),
180            kpp = GET_SLOT(kb, Matrix_pSym);
181      if (ctbl != (nf*(nf + 1))/2)      int *ti = INTEGER(GET_SLOT(tb, Matrix_iSym)),
182          error("length of ctab = %d is not permisable", ctbl);          *tp = INTEGER(tpp),
183      if (j < 1 || j > (nf - 1))          *ii = INTEGER(GET_SLOT(ib, Matrix_iSym)),
184          error("j must be in the range [1,%d]", nf - 1);          *ip = INTEGER(GET_SLOT(ib, Matrix_pSym)),
185      jm1 = j - 1;                /* 0-based index */          *ki = INTEGER(GET_SLOT(kb, Matrix_iSym)),
186      for (i = j; i < nf; i++) {          *kp = INTEGER(kpp),
187          diag_update(ctab, jm1, i);          tnc = length(tpp) - 1,
188          for (k = i; k < nf; k++)          knc = length(kpp) - 1;
189              offdiag_update(ctab, jm1, i, k);      int jj, extra;
190
191        if (k > i || j >= k)
192            error("i,j,k values of %d,%d,%d do not satisfy j < k <= i",
193                  i, j, k);
194                                    /* bound the number of extra elements */
195        extra = 0;
196        for (jj = 0; jj < knc; jj++) {
197            int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1];
198            for (kk = kp[jj]; kk < k2; kk++) {
199                for (i1 = ip[jj]; i1 < i2; i1++) {
200                        if ((check_csc_index(tp, ti, ii[i1], ki[kk]) < 0) &&
201                                    /* only update upper triangle of
202                                     * diagonal blocks */
203                            ((k != i) || (ii[i1] <= ki[kk]))) extra++;
204                }
205            }
206        }
207        if (!extra) return;
208        {
209            int pos, nnz = tp[tnc];
210            int ntot = nnz + extra;
211            int *Ai = Calloc(ntot, int),
212                *Ti = Calloc(ntot, int),
213                *Tj = Calloc(ntot, int),
214                *Dims = INTEGER(GET_SLOT(tb, Matrix_DimSym));
215            double *Ax;
216
217            Memcpy(Ti, ti, nnz);    /* make a copy of the row indices */
218            for (pos = 0, jj = 0; jj < tnc; jj++) { /* fill in the column indices */
219                int j2 = tp[jj + 1];
220                for (; pos < j2; pos++) Tj[pos] = jj;
221            }
222                                    /* add the extra elements */
223            for (jj = 0; jj < knc; jj++) {
224                int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1];
225                for (kk = kp[jj]; kk < k2; kk++) {
226                    for (i1 = ip[jj]; i1 < i2; i1++) {
227                        if ((check_csc_index(tp, ti, ii[i1], ki[kk]) < 0) &&
228                            ((k != i) || (ii[i1] <= ki[kk]))) {
229                            Ti[pos] = ii[i1];
230                            Tj[pos] = ki[kk];
231                            pos++;
232                        }
233                    }
234                }
235            }
236            triplet_to_col(Dims[0], Dims[1], ntot, Ti, Tj, (double *) NULL,
237                           tp, Ai, (double *) NULL);
238            nnz = tp[tnc];
239            SET_SLOT(tb, Matrix_iSym, allocVector(INTSXP, nnz));
240            Memcpy(INTEGER(GET_SLOT(tb, Matrix_iSym)), Ai, nnz);
241            SET_SLOT(tb, Matrix_xSym, allocVector(REALSXP, nnz));
242            Ax = REAL(GET_SLOT(tb, Matrix_xSym));
243            for (j = 0; j < nnz; j++) Ax[j] = 1.;
244            Free(Ai); Free(Ti); Free(Tj);
245            return;
246      }      }

247  }  }
248
249  /**  /**
250   * Permute the levels of one of the grouping factors in a bCrosstab object   * Permute the levels of one of the grouping factors in a bCrosstab object
251   *   *
252   * @param ctab Pointer to a bCrosstab object   * @param ctab Pointer to a bCrosstab object
253   * @param i index (1-based) of the factor levels to permute   * @param nf number of factors in ctab
254     * @param jj index (0-based) of the factor levels to permute
255     * @param ncj number of columns in level jj
256   * @param perm permutation (0-based) to apply   * @param perm permutation (0-based) to apply
257   *   * @param pperm inverse of the permutation
* @return the ctab object with the levels permuted
258   */   */
259  SEXP bCrosstab_permute(SEXP ctab, SEXP ip, SEXP permp)  void bCrosstab_permute(SEXP ctab, int nf, int jj,
260                           const int perm[], const int iperm[])
261  {  {
262      SEXP cs, trip, ipt, jpt;      SEXP trip, ipt, jpt;
263      int *perm, ind = asInteger(ip) - 1,      int j, nnz;
264          ctbl = length(ctab), plen = length(permp), j;
265      int nf = ((int) (-1 + sqrt((double)(1 + 8*ctbl))))/2;      for (j = 0; j < nf; j++) {
266      int *iperm = (int *) Calloc(plen, int);          int ind = (j < jj ? Lind(jj, j) : Lind(j, jj));
267            SEXP csc = VECTOR_ELT(ctab, ind);
268            int *Dims = INTEGER(GET_SLOT(csc, Matrix_DimSym));
269
270            trip = csc_to_triplet(csc);
271            ipt = GET_SLOT(trip, Matrix_iSym);
272            nnz = length(ipt);
273            jpt = GET_SLOT(trip, Matrix_jSym);
274            if (j <= jj) ind_permute(INTEGER(ipt), nnz, iperm);
275            if (j >= jj) ind_permute(INTEGER(jpt), nnz, iperm);
276            if (j == jj)
277                make_upper_triangular(INTEGER(ipt), INTEGER(jpt), nnz);
278            triplet_to_col(Dims[0], Dims[1], nnz, INTEGER(ipt),
279                           INTEGER(jpt),
280                           REAL(GET_SLOT(trip, Matrix_xSym)),
281                           INTEGER(GET_SLOT(csc, Matrix_pSym)),
282                           INTEGER(GET_SLOT(csc, Matrix_iSym)),
283                           REAL(GET_SLOT(csc, Matrix_xSym)));
284        }
285    }
286
287    static void
288    factor_levels_permute(SEXP dest, SEXP src, const int perm[],
289                          const int iperm[])
290    {
291        SEXP dlev = getAttrib(dest, R_LevelsSymbol),
292            slev = getAttrib(src, R_LevelsSymbol);
293        int nlev = length(dlev), flen = length(dest);
294        int *d = INTEGER(dest), *s = INTEGER(src), i;
295
296        if (length(slev) != nlev)
297            error("number of levels in src and dest must match");
298        if (length(src) != flen)
299            error("length of src and dest must match");
300        for (i = 0; i < nlev; i++)
301            SET_STRING_ELT(dlev, i, STRING_ELT(slev, perm[i]));
302        for (i = 0; i < flen; i++)
303            d[i] = 1 + iperm[s[i]-1];
304    }
305
306    SEXP bCrosstab_convert(SEXP bCtab)
307    {
308        char *anms[] = {"flist", "L", "Linv", "perm", "Parent"};
309        SEXP flist = VECTOR_ELT(bCtab, 0),
310            ctab = PROTECT(duplicate(VECTOR_ELT(bCtab, 1))),
311            ans = PROTECT(make_named_list(5, anms));
312        SEXP fcp, perm, L, Linv, Parent;
313        int ctbl = length(ctab), j, nf = length(flist);
314
315      if (ctbl != (nf*(nf + 1))/2)      if (ctbl != (nf*(nf + 1))/2)
316          error("length of ctab = %d is not permisable", ctbl);          error("length of ctab = %d is not permisable", ctbl);
317      if (0 < ind || nf <= ind)      SET_VECTOR_ELT(ans, 0, duplicate(flist));
318          error("i must be in the range [1,%d]", nf);      SET_VECTOR_ELT(ans, 1, allocVector(VECSXP, ctbl));
319      if (!isInteger(permp) || plen !=      SET_VECTOR_ELT(ans, 2, allocVector(VECSXP, nf));
320          (length(GET_SLOT(VECTOR_ELT(ctab, Lind(ind, ind)), Matrix_pSym)) - 1))      SET_VECTOR_ELT(ans, 3, allocVector(VECSXP, nf));
321          error("perm must be an integer vector of length %d",      SET_VECTOR_ELT(ans, 4, allocVector(VECSXP, nf));
322                length(GET_SLOT(VECTOR_ELT(ctab, Lind(ind, ind)), Matrix_pSym)) - 1);      fcp = VECTOR_ELT(ans, 0);
323      if (!ldl_valid_perm(plen, INTEGER(permp), iperm))      L = VECTOR_ELT(ans, 1);
324          error("perm is not a valid 0-based permutation vector");      setAttrib(L, R_NamesSymbol, duplicate(getAttrib(ctab, R_NamesSymbol)));
325      for (j = 0; j < ind; j++) {      Linv = VECTOR_ELT(ans, 2);
326          trip = csc_to_triplet(VECTOR_ELT(ctab, Lind(ind, j)));      setAttrib(Linv, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
327          ipt = GET_SLOT(trip, Matrix_iSym);      perm = VECTOR_ELT(ans, 3);
328          ind_permute(INTEGER(ipt), length(ipt), perm);      setAttrib(perm, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
329          SET_VECTOR_ELT(ctab, Lind(ind, j), triplet_to_csc(trip));      Parent = VECTOR_ELT(ans, 4);
330        setAttrib(Parent, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
331        for (j = 0; j < nf; j++) {
332            int dind = Lind(j, j), i, k;
333            SEXP ctd = VECTOR_ELT(ctab, dind); /* diagonal in crosstab */
334            SEXP Dimslot = GET_SLOT(ctd, Matrix_DimSym),
335                Ljj, cpp = GET_SLOT(ctd, Matrix_pSym),
336                cip = GET_SLOT(ctd, Matrix_iSym);
337            int *Lp, *Perm, *cp = INTEGER(cpp),
338                *ci = INTEGER(cip), *parent,
339                ncj = length(cpp) - 1,
340                nnz = length(cip);
341
342            SET_VECTOR_ELT(Parent, j, allocVector(INTSXP, ncj));
343            parent = INTEGER(VECTOR_ELT(Parent, j));
344            SET_VECTOR_ELT(perm, j, allocVector(INTSXP, ncj));
345            Perm = INTEGER(VECTOR_ELT(perm, j));
346            SET_VECTOR_ELT(L, dind, NEW_OBJECT(MAKE_CLASS("cscMatrix")));
347            Ljj = VECTOR_ELT(L, dind);
348            SET_SLOT(Ljj, Matrix_DimSym, duplicate(Dimslot));
349            SET_SLOT(Ljj, Matrix_factorization, allocVector(VECSXP, 0));
350            SET_SLOT(Ljj, Matrix_pSym, allocVector(INTSXP, ncj + 1));
351            Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));
352            if (nnz > ncj) {        /* calculate fill-reducing permutation */
353                int *iPerm = Calloc(ncj, int),
354                    *Lnz = Calloc(ncj, int);
355                double *Lx;
356
357                ssc_metis_order(ncj, cp, ci, Perm, iPerm);
358                                    /* apply to the crosstabulation */
359                bCrosstab_permute(ctab, nf, j, Perm, iPerm);
360                                    /* apply to the factor */
361                factor_levels_permute(VECTOR_ELT(fcp, j),
362                                      VECTOR_ELT(flist, j), Perm, iPerm);
363                                    /* symbolic analysis to get Parent */
364                ldl_symbolic(ncj, cp, ci, Lp, /* iPerm used as scratch */
365                             parent, Lnz, iPerm,
366                             (int *) NULL, (int *) NULL);
367                nnz = Lp[ncj];
368                Free(iPerm); Free(Lnz);
369                SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));
370                SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, nnz));
371                Lnz = INTEGER(GET_SLOT(Ljj, Matrix_iSym));
372                Lx = REAL(GET_SLOT(Ljj, Matrix_xSym));
373                for (i = 0; i < nnz; i++) {
374                    Lnz[i] = 0; Lx[i] = 0.;
375                }
376            } else {
377                for (i = 0; i <= ncj; i++) {
378                    Lp[i] = 0;
379                    parent[i] = -1;
380                    Perm[i] = i;
381                }
382                SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));
383                SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, 0));
384            }
385            SET_VECTOR_ELT(Linv, j,
386                           Parent_inverse(VECTOR_ELT(Parent, j),
387                                          ScalarLogical(1)));
388    /* FIXME: Update any blocks below the diagonal in this column if necessary*/
389            for (k = j+1; k < nf; k++) { /* update remaining columns, if any */
390                for (i = k; i < nf; i++) block_update(ctab, j, k, i);
391            }
392            for (i= 0; i < j; i++) { /* copy blocks to the left */
393                SET_VECTOR_ELT(L, Lind(j,i), duplicate(VECTOR_ELT(ctab, Lind(j,i))));
394      }      }
trip = csc_to_triplet(VECTOR_ELT(ctab, Lind(ind, ind)));
ipt = GET_SLOT(trip, Matrix_iSym);
ind_permute(INTEGER(ipt), length(ipt), perm);
jpt = GET_SLOT(trip, Matrix_jSym);
ind_permute(INTEGER(jpt), length(jpt), perm);
make_lower_triangular(INTEGER(ipt), INTEGER(jpt), length(jpt));
SET_VECTOR_ELT(ctab, Lind(ind, ind), triplet_to_csc(trip));
for (j = ind + 1; j < nf; j++) {
trip = csc_to_triplet(VECTOR_ELT(ctab, Lind(j, ind)));
jpt = GET_SLOT(trip, Matrix_jSym);
ind_permute(INTEGER(jpt), length(jpt), perm);
SET_VECTOR_ELT(ctab, Lind(j, ind), triplet_to_csc(trip));
395      }      }
396      return ctab;      UNPROTECT(2);
397        return ans;
398  }  }
399
400
401
402

Legend:
 Removed from v.347 changed lines Added in v.356