SCM

SCM Repository

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

Diff of /pkg/src/ssclme.c

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

revision 125, Sun Apr 25 11:59:37 2004 UTC revision 245, Mon Jul 12 12:46:54 2004 UTC
# Line 1  Line 1 
1  #include "ssclme.h"  #include "ssclme.h"
2    
3  /**  #define slot_dup(dest, src, sym)  SET_SLOT(dest, sym, duplicate(GET_SLOT(src, sym)))
  * 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;  
 }  
4    
5  /**  /**
6   * Determine if a fill-reducing permutation is needed for the pairwise   * Using the sscCrosstab object from the grouping factors, generate
7   * cross-tabulation matrix.  If so, determine such a permutation   * the slots in an ssclme object related to the symmetric sparse
8   * (using Metis) then separate the groups.   * matrix representation of Z'Z.  If the model matrices for the
9   *   * grouping factors have only one column each then the structure can
10   * @param ctab pointer to a pairwise cross-tabulation object   * be copied, otherwise it must be generated from the sscCrosstab and
11   *   * the number of columns per grouping factor.
12   * @return pointer to an integer R vector.   *
13     * @param nf number of factors
14     * @param nc vector of length nf+2 with number of columns in model matrices
15     * @param ctab pointer to the sscCrosstab object
16     * @param ssc pointer to an ssclme object to be filled out
17   */   */
   
 SEXP ctab_permute(SEXP ctab)  
 {  
     SEXP val, GpSl = GET_SLOT(ctab, Matrix_GpSym);  
     int *Ai = INTEGER(GET_SLOT(ctab, Matrix_iSym)),  
         *Ap = INTEGER(GET_SLOT(ctab, Matrix_pSym)),  
         *Gp = INTEGER(GpSl),  
         *perm,  
         *work,  
         i,  
         j,  
         n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],  
         nf = length(GpSl) - 1,  
         pos;  
   
     if (ctab_isNested(n, nf, 1, Ap, Ai, Gp))  
         return allocVector(INTSXP, 0);  
     val =  allocVector(INTSXP, n);  
     perm = INTEGER(val);  
     work = (int *) R_alloc(n, sizeof(int));  
     ssc_metis_order(n, Ap, Ai, work, perm);     /* perm gets inverse perm */  
     /* work now contains desired permutation but with groups scrambled */  
   
     /* copy work into perm preserving the order of the groups */  
     pos = 0;            /* position in new permutation */  
     for (i = 0; i < nf; i++) {  
         for (j = 0; j < n; j++) {  
             int jj = work[j];  
             if (Gp[i] <= jj && jj < Gp[i+1]) {  
                 perm[pos] = jj;  
                 pos++;  
             }  
         }  
     }  
     return val;  
 }  
   
18  static  static
19  void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)  void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)
20  {  {
21      int *snc, i, copyonly = 1;      int *snc, i, copyonly = 1;
22    
23      for (i = 0; i < nf; i++) {      SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));
24          if (nc[i] > 1) copyonly = 0;      snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
25        for (i = 0; i <= nf; i++) {
26            snc[i] = nc[i];
27            if (nc[i] > 1 && i < nf) copyonly = 0;
28      }      }
29      if (copyonly) {      if (copyonly) {
30          SET_SLOT(ssc, Matrix_pSym, duplicate(GET_SLOT(ctab, Matrix_pSym)));          slot_dup(ssc, ctab, Matrix_pSym);
31          SET_SLOT(ssc, Matrix_iSym, duplicate(GET_SLOT(ctab, Matrix_iSym)));          slot_dup(ssc, ctab, Matrix_iSym);
32          SET_SLOT(ssc, Matrix_xSym, duplicate(GET_SLOT(ctab, Matrix_xSym)));          slot_dup(ssc, ctab, Matrix_xSym);
33          SET_SLOT(ssc, Matrix_DimSym,          slot_dup(ssc, ctab, Matrix_DimSym);
34                   duplicate(GET_SLOT(ctab, Matrix_DimSym)));          slot_dup(ssc, ctab, Matrix_GpSym);
35          SET_SLOT(ssc, Matrix_GpSym, duplicate(GET_SLOT(ctab, Matrix_GpSym)));          return;
36      } else {      }
37        {
38          int          int
39              *GpIn = INTEGER(GET_SLOT(ctab, Matrix_GpSym)),              *GpIn = INTEGER(GET_SLOT(ctab, Matrix_GpSym)),
40              *GpOut,              *GpOut,
# Line 130  Line 60 
60              }              }
61          }          }
62          nOut = GpOut[nf];       /* size of output matrix */          nOut = GpOut[nf];       /* size of output matrix */
63          SET_SLOT(ssc, Matrix_DimSym,          SET_SLOT(ssc, Matrix_DimSym, allocVector(INTSXP, 2));
                  duplicate(GET_SLOT(ctab, Matrix_DimSym)));  
64          dims = INTEGER(GET_SLOT(ssc, Matrix_DimSym));          dims = INTEGER(GET_SLOT(ssc, Matrix_DimSym));
65          dims[0] = dims[1] = nOut;          dims[0] = dims[1] = nOut;
66          SET_SLOT(ssc, Matrix_pSym, allocVector(INTSXP, nOut + 1));          SET_SLOT(ssc, Matrix_pSym, allocVector(INTSXP, nOut + 1));
# Line 162  Line 91 
91          for (i = 0; i < nf; i++) { /* fill in the rows */          for (i = 0; i < nf; i++) { /* fill in the rows */
92              int j, jj, nci = nc[i], p2 = GpIn[i+1];              int j, jj, nci = nc[i], p2 = GpIn[i+1];
93              for (j = GpIn[i]; j < p2; j++) { /* first col for input col */              for (j = GpIn[i]; j < p2; j++) { /* first col for input col */
94                  int ii = AiIn[j], mj = map[j], ncci = ncc[ii],                  int k, mj = map[j], p3 = ApIn[j+1], pos = ApOut[mj];
95                      pos = ApOut[mj];                  for (k = ApIn[j]; k < p3; k++) {
96                        int ii = AiIn[k], ncci = ncc[ii];
97                  AiOut[pos++] = map[ii];                  AiOut[pos++] = map[ii];
98                  if (ii < j) {   /* above the diagonal */                  if (ii < j) {   /* above the diagonal */
99                      for (jj = 1; jj < ncci; jj++) {                      for (jj = 1; jj < ncci; jj++) {
# Line 175  Line 105 
105                                                  * another diagonal element */                                                  * another diagonal element */
106                      int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];                      int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];
107                      Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);                      Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);
108                      AiOut[ApOut[mjj + 1] - 1] = mjj; /* maybe mjj-1? */                          AiOut[ApOut[mjj + 1] - 1] = mjj;
109                  }                  }
110              }              }
111          }          }
         Free(map); Free(ncc);  
112      }      }
113      SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));          Free(map); Free(ncc);
     snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));  
     for (i = 0; i <= nf; i++) {  
         snc[i] = nc[i];  
114      }      }
115  }  }
116    
117  void ssclme_fill_LIp(int n, const int Parent[], int LIp[])  /**
118     * Calculate and store the maximum number of off-diagonal elements in
119     * the inverse of L, based on the elimination tree.  The maximum is
120     * itself stored in the Parent array.  (FIXME: come up with a better design.)
121     *
122     * @param n number of columns in the matrix
123     * @param Parent elimination tree for the matrix
124     */
125    static void ssclme_calc_maxod(int n, int Parent[])
126  {  {
127      int *sz = Calloc(n, int), i;      int *sz = Calloc(n, int), i, mm = -1;
128      for (i = n - 1; i >= 0; i--) {      for (i = n - 1; i >= 0; i--) {
129          sz[i] = (Parent[i] < 0) ? 0 : 1 + sz[Parent[i]];          sz[i] = (Parent[i] < 0) ? 0 : (1 + sz[Parent[i]]);
130            if (sz[i] > mm) mm = sz[i];
131      }      }
132      LIp[0] = 0;      Parent[n] = mm;
     for (i = 0; i < n; i++) LIp[i+1] = LIp[i] + sz[i];  
133      Free(sz);      Free(sz);
134  }  }
135    
136  static  /**
137  void ssclme_fill_LIi(int n, const int Parent[], const int LIp[], int LIi[])   * Create an ssclme object from a list of grouping factors, sorted in
138  {   * order of non-increasing numbers of levels, and an integer vector of
139      int i;   * the number of columns in the model matrices.  There is one more
140      for (i = n; i > 0; i--) {   * element in ncv than in facs.  The last element is the number of
141          int im1 = i - 1, Par = Parent[im1];   * columns in the model matrix for the fixed effects plus the
142          if (Par >= 0) {   * response.  (i.e. p+1)
143              LIi[LIp[im1]] = Par;   *
144              Memcpy(LIi + LIp[im1] + 1, LIi + LIp[Par],   * @param facs pointer to a list of grouping factors
145                     LIp[Par + 1] - LIp[Par]);   * @param ncv pointer to an integer vector of number of columns per model matrix
146          }   *
147      }   * @return pointer to an ssclme object
148  }   */
   
149  SEXP  SEXP
150  ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)  ssclme_create(SEXP facs, SEXP ncv)
151  {  {
152      SEXP ctab, nms, ssc, tmp,      SEXP ctab, nms, ssc, tmp,
153          val = PROTECT(allocVector(VECSXP, 2)),          val = PROTECT(allocVector(VECSXP, 2)),
154          dd = PROTECT(allocVector(INTSXP, 3));   /* dimensions of 3-D arrays */          dd = PROTECT(allocVector(INTSXP, 3));   /* dimensions of 3-D arrays */
155      int *Ai, *Ap, *Gp, *LIp, *Lp, *Parent,      int *Ai, *Ap, *Gp, *Lp, *Parent,
156          *nc, Lnz, i, nf = length(facs), nzcol, pp1,          *nc, Lnz, i, nf = length(facs), nzcol, pp1,
157          *dims = INTEGER(dd);          *dims = INTEGER(dd);
158    
# Line 230  Line 163 
163      ssc = VECTOR_ELT(val, 0);      ssc = VECTOR_ELT(val, 0);
164                                  /* Pairwise cross-tabulation */                                  /* Pairwise cross-tabulation */
165      ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));      ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));
166      SET_VECTOR_ELT(val, 1, ctab_permute(ctab));      SET_VECTOR_ELT(val, 1, sscCrosstab_groupedPerm(ctab));
167      if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */      if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */
168          ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],          ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
169                               1, INTEGER(VECTOR_ELT(val, 1)),                               1, INTEGER(VECTOR_ELT(val, 1)),
# Line 259  Line 192 
192             sizeof(double) * pp1 * pp1);             sizeof(double) * pp1 * pp1);
193      SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));      SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));
194      Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));      Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));
195      SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol));      SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol + 1));
196      Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));      Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));
197      SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));      SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));
198      SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));      SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));
# Line 267  Line 200 
200                   (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */                   (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */
201                   (int *) R_alloc(nzcol, sizeof(int)), /* Flag */                   (int *) R_alloc(nzcol, sizeof(int)), /* Flag */
202                   (int *) NULL, (int *) NULL); /* P & Pinv */                   (int *) NULL, (int *) NULL); /* P & Pinv */
203        ssclme_calc_maxod(nzcol, Parent);
204      Lnz = Lp[nzcol];      Lnz = Lp[nzcol];
205      SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));      SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));
206      SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));      SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));
# Line 304  Line 238 
238          memset(REAL(VECTOR_ELT(tmp, i)), 0,          memset(REAL(VECTOR_ELT(tmp, i)), 0,
239                 sizeof(double) * nci * nci * mi);                 sizeof(double) * nci * nci * mi);
240      }      }
     SET_SLOT(ssc, Matrix_LIpSym, allocVector(INTSXP, nzcol + 1));  
     LIp = INTEGER(GET_SLOT(ssc, Matrix_LIpSym));  
     ssclme_fill_LIp(nzcol, Parent, LIp);  
     if (asInteger(threshold) > (Lnz = LIp[nzcol])) {  
         SET_SLOT(ssc, Matrix_LIiSym, allocVector(INTSXP, Lnz));  
         ssclme_fill_LIi(nzcol, Parent, LIp,  
                         INTEGER(GET_SLOT(ssc, Matrix_LIiSym)));  
         SET_SLOT(ssc, Matrix_LIxSym, allocVector(REALSXP, Lnz));  
         memset(REAL(GET_SLOT(ssc, Matrix_LIxSym)), 0,  
                sizeof(double) * Lnz);  
     }  
241      UNPROTECT(2);      UNPROTECT(2);
242      return val;      return val;
243  }  }
244    
245    /**
246     * Copy information on Z'Z accumulated in the bVar array to Z'Z
247     *
248     * @param ncj number of columns in this block
249     * @param Gpj initial column for this group
250     * @param Gpjp initial column for the next group
251     * @param bVj pointer to the ncj x ncj x mj array to be filled
252     * @param Ap column pointer array for Z'Z
253     * @param Ai row indices for Z'Z
254     * @param Ax elements of Z'Z
255     */
256  static  static
257  void bVj_to_A(int ncj, int Gpj, int Gpjp, const double bVj[],  void bVj_to_A(int ncj, int Gpj, int Gpjp, const double bVj[],
258                const int Ap[], const int Ai[], double Ax[])                const int Ap[], const int Ai[], double Ax[])
# Line 335  Line 269 
269      }      }
270  }  }
271    
272    /**
273     * Copy the dimnames from the list of grouping factors and the model
274     * matrices for the grouping factors into the appropriate parts of the
275     * ssclme object.
276     *
277     * @param x pointer to an ssclme object
278     * @param facs pointer to a list of factors
279     * @param mmats pointer to a list of model matrices
280     *
281     * @return NULL
282     */
283  SEXP  SEXP
284  ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)  ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)
285  {  {
# Line 360  Line 305 
305      return R_NilValue;      return R_NilValue;
306  }  }
307    
308    /**
309     * Update the numerical entries x, ZtX, and XtX in an ssclme object
310     * according to a set of model matrices.
311     *
312     * @param x pointer to an ssclme object
313     * @param facs pointer to a list of grouping factors
314     * @param mmats pointer to a list of model matrices
315     *
316     * @return NULL
317     */
318  SEXP  SEXP
319  ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)  ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)
320  {  {
# Line 369  Line 324 
324          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
325          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
326          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
327            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
328          i, j, k,          i, j, k,
329          ione = 1,          ione = 1,
330          nf = length(mmats) - 1,          nf = length(mmats) - 1,
# Line 435  Line 391 
391                  nck = nc[k],                  nck = nc[k],
392                  scalar = ncj == 1 && nck == 1;                  scalar = ncj == 1 && nck == 1;
393              double              double
394                  *Zk = Z[k], *work;                  *Zk = Z[k], *work = (double *) NULL;
395    
396              if (!scalar) work = Calloc(ncj * nck, double);              if (!scalar) work = Calloc(ncj * nck, double);
397              for (i = 0; i < nobs; i++) {              for (i = 0; i < nobs; i++) {
398                  int ii, ind = -1, fpji = fpj[i] - 1,                  int ii, ind = -1, fpji = fpj[i] - 1,
# Line 471  Line 428 
428      }      }
429      Free(Z);      Free(Z);
430      ssclme_transfer_dimnames(x, facs, mmats);      ssclme_transfer_dimnames(x, facs, mmats);
431        status[0] = status[1] = 0;
432      return R_NilValue;      return R_NilValue;
433  }  }
434    
435    /**
436     * Inflate Z'Z according to Omega and create the factorization LDL'
437     *
438     * @param x pointer to an ssclme object
439     *
440     * @return NULL
441     */
442  SEXP ssclme_inflate_and_factor(SEXP x)  SEXP ssclme_inflate_and_factor(SEXP x)
443  {  {
444      SEXP      SEXP
# Line 528  Line 493 
493      return R_NilValue;      return R_NilValue;
494  }  }
495    
496    
497    /**
498     * If status[["factored"]] is FALSE, create and factor Z'Z+Omega, then
499     * create RZX and RXX, the deviance components, and the value of the
500     * deviance for both ML and REML.
501     *
502     * @param x pointer to an ssclme object
503     *
504     * @return NULL
505     */
506  SEXP ssclme_factor(SEXP x)  SEXP ssclme_factor(SEXP x)
507  {  {
508      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
# Line 621  Line 596 
596      return R_NilValue;      return R_NilValue;
597  }  }
598    
599    /**
600     * Return the position of probe in the sorted index vector ind.  It is
601     * known that the position is greater than or equal to start so a linear
602     * search from start is used.
603     *
604     * @param probe value to be matched
605     * @param start index at which to start
606     * @param ind vector of indices
607     *
608     * @return index of the entry matching probe
609     */
610  static  static
611  int ldl_update_ind(int probe, int start, const int ind[])  int ldl_update_ind(int probe, int start, const int ind[])
612  {  {
# Line 630  Line 616 
616  }  }
617    
618  /**  /**
619   * Create the inverse of L and update the diagonal blocks of the inverse   * Update the diagonal blocks of the inverse of LDL' (=Z'Z+W).  The
620   * of LDL' (=Z'Z+W)   * lower Cholesky factors of the updated blocks are stored in the bVar
621     * slot.
622   *   *
623   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
624   *   *
# Line 639  Line 626 
626    
627   */   */
628  static  static
629  SEXP ldl_inverse(SEXP x)  void ldl_inverse(SEXP x)
630  {  {
631      SEXP      SEXP
632          Gpsl = GET_SLOT(x, Matrix_GpSym),          Gpsl = GET_SLOT(x, Matrix_GpSym),
         LIisl = GET_SLOT(x, Matrix_LIiSym),  
         LIpsl = GET_SLOT(x, Matrix_LIpSym),  
633          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
634      int *Gp = INTEGER(Gpsl),      int *Gp = INTEGER(Gpsl),
635          *Li,          *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym)),
636          *LIp = INTEGER(LIpsl), *Lp,          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
637          i,          i,
638          nf = length(Gpsl) - 1,          nf = length(Gpsl) - 1,
639          nzc = length(LIpsl) - 1;          nzc = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
640      double      int maxod = Parent[nzc];
641          *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),      double *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym));
         *Lx;  
642    
643      ssclme_factor(x);      ssclme_factor(x);
644      if (LIp[nzc] == 0) {        /* L and LI are the identity */      if (maxod == 0) {           /* L and L^{-1} are the identity */
645          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
646              Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],              Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],
647                     Gp[i+1] - Gp[i]);                     Gp[i+1] - Gp[i]);
648          }          }
649          return R_NilValue;      } else {
650      }          int *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
651      Lp = INTEGER(GET_SLOT(x, Matrix_LpSym));              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym));
     Li = INTEGER(GET_SLOT(x, Matrix_LiSym));  
     Lx = REAL(GET_SLOT(x, Matrix_LxSym));  
     if (length(LIisl) == LIp[nzc]) { /* LIi is filled */  
         int *LIi = INTEGER(LIisl),  
             *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),  
             j, jj, k, kk, p1, p2, pi1, pi2;  
   
         double *LIx = REAL(GET_SLOT(x, Matrix_LIxSym)),  
             one = 1., zero = 0.;  
   
         memset(LIx, 0, sizeof(double) * LIp[nzc]);  
                                 /* calculate inverse */  
         for (i = 0; i < nzc; i++) {  
             p1 = Lp[i]; p2 = Lp[i+1]; pi1 = LIp[i]; pi2 = LIp[i+1];  
                                 /* initialize from unit diagonal term */  
             kk = pi1;  
             for (j = p1; j < p2; j++) {  
                 k = Li[j];  
                 while (LIi[kk] < k && kk < pi2) kk++;  
                 if (LIi[kk] != k) error("logic error in ldl_inverse");  
                 LIx[kk] = -Lx[j];  
             }  
             for (j = pi1; j < pi2; j++) {  
                 jj = LIi[j];  
                 p1 = Lp[jj]; p2 = Lp[jj+1];  
                 kk = j;  
                 for (jj = p1; jj < p2; jj++) {  
                     k = Li[jj];  
                     while (LIi[kk] < k && kk < pi2) kk++;  
                     if (LIi[kk] != k) error("logic error in ldl_inverse");  
                     LIx[kk] -= Lx[jj]*LIx[j];  
                 }  
             }  
         }  
         for (i = 0; i < nf; i++) { /* accumulate bVar */  
             int G1 = Gp[i], G2 = Gp[i+1], j, k, kk,  
                 nci = nc[i], nr, nr1, rr;  
             double *bVi = REAL(VECTOR_ELT(bVar, i)), *tmp;  
652    
653              nr = -1;          double one = 1.0, zero = 0.,
654              for (j = G1; j < G2; j += nci) {              *Lx = REAL(GET_SLOT(x, Matrix_LxSym));
                 rr = 1 + LIp[j + 1] - LIp[j];  
                 if (rr > nr) nr = rr;  
             }  
             tmp = Calloc(nr * nci, double); /* scratch storage */  
             nr1 = nr + 1;  
                                 /* initialize bVi to zero */  
             memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);  
             for (j = G1; j < G2; j += nci) {  
                 memset(tmp, 0, sizeof(double) * nr * nci);  
                 rr = 1 + LIp[j + 1] - LIp[j];  
                 for (k = 0; k < nci; k++) { /* copy columns */  
                     tmp[k * nr1] = 1.; /* (unstored) diagonal elt  */  
                     Memcpy(tmp + k*nr1 + 1, LIx + LIp[j + k], rr - k - 1);  
                 }  
                                 /* scale the rows */  
                 tmp[0] = DIsqrt[j]; /* first row only has one non-zero */  
                 for (kk = 1; kk < rr; kk++) {  
                     for (k = 0; k < nci; k++) {  
                         tmp[k * nr + kk] *= DIsqrt[LIi[LIp[j] + kk - 1]];  
                     }  
                 }  
                 F77_CALL(dsyrk)("L", "T", &nci, &rr, &one, tmp, &nr,  
                                 &zero, bVi + (j - G1) * nci, &nci);  
                 F77_CALL(dpotrf)("L", &nci, bVi + (j - G1) * nci,  
                                  &nci, &kk);  
                 if (kk)         /* should never happen */  
                     error(  
                         "Rank deficient variance matrix at group %d, level %d",  
                         i + 1, j + 1);  
             }  
             Free(tmp);  
         }  
         return R_NilValue;  
     }  
     if (length(LIisl)) error("logic error in ssclme_ldl_inverse");  
     else {                  /* LIi and LIx are too big and not used */  
         int *counts = Calloc(nzc, int), info, maxod = -1;  
         int *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym));  
         int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));  
         double one = 1.0, zero = 0.;  
                                 /* determine maximum # of off-diagonals */  
         for (i = nzc - 1; i >= 0; i--) { /* in a column of L^{-1} */  
             counts[i] = (Parent[i] < 0) ? 0 : 1 + counts[Parent[i]];  
             if (counts[i] > maxod) maxod = counts[i];  
         }  
         Free(counts);  
655    
656          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
657              int j, jj, k, kk, nci = nc[i], nr, p, p2, pp,              int j, jj, k, kk, nci = nc[i], nr, p, p2, pj, pp,
658                  m = maxod + nci,                  m = maxod + 1,
659                  *ind = Calloc(m, int);                  *ind = Calloc(m, int), G1 = Gp[i], G2 = Gp[i+1];
660              double              double
661                  *tmp = Calloc(m * nci, double),                  *tmp = Calloc(m * nci, double),
662                  *mpt = REAL(VECTOR_ELT(bVar, i));                  *bVi = REAL(VECTOR_ELT(bVar, i));
663    
664              for (j = Gp[i]; j < Gp[i+1]; j += nci) {                                  /* initialize bVi to zero */
665                  memset((void *) tmp, 0, sizeof(double) * m * nci);              memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
666    
667                  kk = 0;         /* ind holds indices of non-zeros */              for (j = G1; j < G2; j += nci) {
668                    kk = 0;         /* ind gets indices of non-zeros */
669                  jj = j;         /* in this block of columns */                  jj = j;         /* in this block of columns */
670                  while (jj >= 0) {                  while (jj >= 0) {
671                      ind[kk++] = jj;                      ind[kk++] = jj;
# Line 776  Line 677 
677                  for (k = 0; k < nci; k++) {                  for (k = 0; k < nci; k++) {
678                      double *ccol = tmp + k * nr;                      double *ccol = tmp + k * nr;
679    
680                      ccol[k] = 1.;                      for (kk = 0; kk < nr; kk++) ccol[kk] = 0.;
681                      kk = k;                      ccol[k] = 1.; /* initialize from unit diagonal */
682                      for (jj = j + k; jj >= 0; jj = Parent[jj]) {                      for (jj = j + k; jj >= 0; jj = Parent[jj]) {
683                          p2 = Lp[jj+1];                          p2 = Lp[jj+1];
684                          pp = kk;                          pp = pj = ldl_update_ind(jj, 0, ind);
685                          for (p = Lp[jj]; p < p2; p++) {                          for (p = Lp[jj]; p < p2; p++) {
686                              pp = ldl_update_ind(Li[p], pp, ind);                              pp = ldl_update_ind(Li[p], pp, ind);
687                              ccol[pp] -= Lx[p] * ccol[kk];                              ccol[pp] -= Lx[p] * ccol[pj];
688                          }                          }
689                      }                      }
690                  }                  }
691                                  /* scale rows */  
692                  for (kk = 0; kk < nr; kk++) {                  for (kk = 0; kk < nr; kk++) { /* scale rows */
693                      for (k = 0; k < nci; k++) {                      for (k = 0; k < nci; k++) {
694                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];
695                      }                      }
696                  }                  }
697                  F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,                  F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,
698                                  &zero, mpt + (j - Gp[i])*nci, &nci);                                  &zero, bVi + (j - G1)*nci, &nci);
699                  F77_CALL(dpotrf)("L", &nci, mpt + (j - Gp[i])*nci,                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,
700                                   &nci, &info);                                   &nci, &jj);
701                  if (info)       /* should never happen */                  if (jj)         /* should never happen */
702                      error(                      error(
703                          "Rank deficient variance matrix at group %d, level %d",                          "Rank deficient variance matrix at group %d, level %d, error code %d",
704                          i + 1, j + 1);                          i + 1, j + 1, jj);
705              }              }
706              Free(tmp); Free(ind);              Free(tmp); Free(ind);
707          }          }
708      }      }
     return R_NilValue;  
709  }  }
710    
711    /**
712     * If necessary, factor Z'Z+Omega, ZtX, and XtX then, if necessary,
713     * form RZX, RXX, and bVar for the inverse of the Cholesky factor.
714     *
715     * @param x pointer to an ssclme object
716     *
717     * @return NULL (x is updated in place)
718     */
719  SEXP ssclme_invert(SEXP x)  SEXP ssclme_invert(SEXP x)
720  {  {
721      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
# Line 847  Line 755 
755      return R_NilValue;      return R_NilValue;
756  }  }
757    
758    /**
759     * Create and insert initial values for Omega_i.
760     *
761     * @param x pointer to an ssclme object
762     *
763     * @return NULL
764     */
765  SEXP ssclme_initial(SEXP x)  SEXP ssclme_initial(SEXP x)
766  {  {
767      SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),      SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),
# Line 971  Line 886 
886                                          nobs + 1 - pp1 : nobs))));                                          nobs + 1 - pp1 : nobs))));
887  }  }
888    
889    /**
890     * Calculate the length of the parameter vector, which is called coef
891     * for historical reasons.
892     *
893     * @param nf number of factors
894     * @param nc number of columns in the model matrices for each factor
895     *
896     * @return total length of the coefficient vector
897     */
898  static  static
899  int coef_length(int nf, const int nc[])  int coef_length(int nf, const int nc[])
900  {  {
# Line 980  Line 904 
904  }  }
905    
906  /**  /**
907   * Extract the upper triangles of the Omega matrices.   * Extract the upper triangles of the Omega matrices.  These aren't
908   * (These aren't "coefficients" but the extractor is   * "coefficients" but the extractor is called coef for historical
909   * called coef for historical reasons.)   * reasons.  Within each group these values are in the order of the
910     * diagonal entries first then the strict upper triangle in row
911     * order.
912   *   *
913   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
914   *   *
# Line 1021  Line 947 
947    
948  /**  /**
949   * Extract the unconstrained parameters that determine the   * Extract the unconstrained parameters that determine the
950   * Omega matrices. (Called coef for historical reasons.)   * Omega matrices. (Called coef for historical reasons.)  The
951     * unconstrained parameters are derived from the LDL' decomposition of
952     * Omega_i.  The first nc[i] entries in each group are the diagonals
953     * of log(D) followed by the strict lower triangle of L in column
954     * order.
955   *   *
956   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
957   *   *
# Line 1161  Line 1091 
1091      return x;      return x;
1092  }  }
1093    
1094  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)  
1095    /**
1096     * Returns the inverse of the updated Omega matrices for an ECME
1097     * iteration.  These matrices are also used in the gradient calculation.
1098     *
1099     * @param x pointer to an ssclme object
1100     * @param REML indicator of REML being used
1101     * @param val pointer to a list of symmetric q_i by q_i matrices
1102     */
1103    static void
1104    common_ECME_gradient(SEXP x, int REML, SEXP val)
1105  {  {
1106      SEXP      SEXP
         Omega = GET_SLOT(x, Matrix_OmegaSym),  
1107          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ncsl = GET_SLOT(x, Matrix_ncSym),  
1108          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1109      int      int
1110          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1111          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1112          *nc = INTEGER(ncsl),          i, j, n = *INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1113          *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),          nf = length(val), nobs = nc[nf + 1], p = nc[nf] - 1;
         REML = asLogical(REMLp),  
         i, info, iter,  
         n = dims[0],  
         nEM = asInteger(nsteps),  
         nf = length(ncsl) - 2,  
         nobs = nc[nf + 1],  
         p,  
         pp1 = dims[1],  
         verbose = asLogical(verb);  
1114      double      double
1115          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
1116          *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),          *b = RZX + p * INTEGER(getAttrib(RZXsl, R_DimSymbol))[0],
1117          *b,          one = 1.0, zero = 0.0;
         alpha,  
         one = 1.,  
         zero = 0.;  
1118    
     p = pp1 - 1;  
     b = RZX + p * n;  
     if (verbose) Rprintf("  EM iterations\n");  
     for (iter = 0; iter <= nEM; iter++) {  
1119          ssclme_invert(x);          ssclme_invert(x);
         if (verbose) {  
             SEXP coef = PROTECT(ssclme_coef(x));  
             int lc = length(coef); double *cc = REAL(coef);  
             Rprintf("%3d %.3f", iter, dev[REML ? 1 : 0]);  
             for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);  
             Rprintf("\n");  
             UNPROTECT(1);  
         }  
1120          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
1121              int ki = Gp[i+1] - Gp[i],              int ki = Gp[i+1] - Gp[i],
1122                  nci = nc[i],                  nci = nc[i],
1123                  mi = ki/nci;                  mi = ki/nci;
1124              double              double
1125                  *vali = REAL(VECTOR_ELT(Omega, i));              *vali = REAL(VECTOR_ELT(val, i)),
1126                alpha = (double)(REML ? (nobs-p) : nobs)/((double) mi);
1127    
             alpha = ((double)(REML?(nobs-p):nobs))/((double)mi);  
1128              F77_CALL(dsyrk)("U", "N", &nci, &mi,              F77_CALL(dsyrk)("U", "N", &nci, &mi,
1129                              &alpha, b + Gp[i], &nci,                              &alpha, b + Gp[i], &nci,
1130                              &zero, vali, &nci);                              &zero, vali, &nci);
# Line 1219  Line 1133 
1133                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1134                              &one, vali, &nci);                              &one, vali, &nci);
1135              if (REML) {              if (REML) {
                 int j;  
1136                  for (j = 0; j < p; j++) {                  for (j = 0; j < p; j++) {
1137                      F77_CALL(dsyrk)("U", "N", &nci, &mi,                      F77_CALL(dsyrk)("U", "N", &nci, &mi,
1138                                  &alpha, RZX + Gp[i] + j*n, &nci,                                  &alpha, RZX + Gp[i] + j*n, &nci,
1139                                  &one, vali, &nci);                                  &one, vali, &nci);
1140                  }                  }
1141              }              }
1142              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);      }
1143    }
1144    
1145    /**
1146     * Perform a number of ECME steps for the REML or ML criterion.
1147     *
1148     * @param x pointer to an ssclme object
1149     * @param nsteps pointer to an integer scalar giving the number of ECME steps to perform
1150     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1151     * @param verb pointer to a logical scalar indicating verbose mode
1152     *
1153     * @return NULL
1154     */
1155    SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1156    {
1157        SEXP
1158            Omega = GET_SLOT(x, Matrix_OmegaSym);
1159        int
1160            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1161            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1162            REML = asLogical(REMLp),
1163            i, info, iter,
1164            nEM = asInteger(nsteps),
1165            nf = length(Omega),
1166            verbose = asLogical(verb);
1167        double
1168            *dev = REAL(GET_SLOT(x, Matrix_devianceSym));
1169    
1170        if (verbose) {
1171            SEXP coef = PROTECT(ssclme_coef(x));
1172            int lc = length(coef); double *cc = REAL(coef);
1173    
1174            ssclme_factor(x);
1175            Rprintf("  EM iterations\n");
1176            Rprintf("%3d %.3f", 0, dev[REML ? 1 : 0]);
1177            for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1178            Rprintf("\n");
1179            UNPROTECT(1);
1180        }
1181        for (iter = 0; iter < nEM; iter++) {
1182            common_ECME_gradient(x, REML, Omega);
1183            status[0] = status[1] = 0;
1184            for (i = 0; i < nf; i++) {
1185                double *vali = REAL(VECTOR_ELT(Omega, i));
1186    
1187                F77_CALL(dpotrf)("U", nc + i, vali, nc + i, &info);
1188              if (info)              if (info)
1189                  error("DPOTRF returned error code %d in Omega[[%d]] update",                  error("DPOTRF returned error code %d in Omega[[%d]] update",
1190                        info, i + 1);                        info, i + 1);
1191              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);              F77_CALL(dpotri)("U", nc + i, vali, nc + i, &info);
1192              if (info)              if (info)
1193                  error("DPOTRI returned error code %d in Omega[[%d]] update",                  error("DPOTRI returned error code %d in Omega[[%d]] update",
1194                        info, i + 1);                        info, i + 1);
1195          }          }
1196          status[0] = status[1] = 0;          if (verbose) {
1197                SEXP coef = PROTECT(ssclme_coef(x));
1198                int lc = length(coef); double *cc = REAL(coef);
1199    
1200                ssclme_factor(x);
1201                Rprintf("%3d %.3f", iter + 1, dev[REML ? 1 : 0]);
1202                for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1203                Rprintf("\n");
1204                UNPROTECT(1);
1205            }
1206      }      }
1207      ssclme_factor(x);      ssclme_factor(x);
1208      return R_NilValue;      return R_NilValue;
1209  }  }
1210    
1211    /**
1212     * Evaluate the gradient with respect to the indicators of the
1213     * positions in the Omega matrices.
1214     *
1215     * @param x Pointer to an ssclme object
1216     * @param REML indicator of whether REML is to be used
1217     * @param val Pointer to an object of the same structure as Omega
1218     */
1219    static void indicator_gradient(SEXP x, int REML, SEXP val)
1220    {
1221        SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1222        int
1223            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1224            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1225            i, nf = length(Omega);
1226    
1227        common_ECME_gradient(x, REML, val);
1228        for (i = 0; i < nf; i++) {
1229            int info, nci = nc[i], mi = (Gp[i+1] - Gp[i])/nci;
1230            double
1231                *work = Memcpy((double *) Calloc(nci * nci, double),
1232                               REAL(VECTOR_ELT(Omega, i)), nci * nci),
1233                alpha = (double) -mi, beta = (double) mi;
1234    
1235            F77_CALL(dpotrf)("U", &nci, work, &nci, &info);
1236            if (info)
1237                error("DPOTRF gave error code %d on Omega[[%d]]", info, i + 1);
1238            F77_CALL(dtrtri)("U", "N", &nci, work, &nci, &info);
1239            if (info)
1240                error("DPOTRI gave error code %d on Omega[[%d]]", info, i + 1);
1241            F77_CALL(dsyrk)("U", "N", &nci, &nci, &alpha, work, &nci,
1242                            &beta, REAL(VECTOR_ELT(val, i)), &nci);
1243            Free(work);
1244        }
1245    }
1246    
1247    /**
1248     * Overwrite a gradient with respect to positions in Omega[[i]] by the
1249     * gradient with respect to the unconstrained parameters.
1250     *
1251     * @param grad pointer to a gradient wrt positions.  Overwritten.
1252     * @param Omega pointer to a list of symmetric (upper storage)
1253     * matrices Omega[[i]].
1254     */
1255    static void unconstrained_gradient(SEXP grad, SEXP Omega)
1256    {
1257        int i, ii, j, nf = length(Omega);
1258        double one = 1.0, zero = 0.0;
1259    
1260        for (i = 0; i < nf; i++) {
1261            SEXP Omegai = VECTOR_ELT(Omega, i);
1262            int nci = *INTEGER(getAttrib(Omegai, R_DimSymbol)),
1263                ncisq = nci * nci, ncip1 = nci + 1;
1264            double *chol =
1265                Memcpy(Calloc(ncisq, double), REAL(Omegai), ncisq),
1266                *ai = REAL(VECTOR_ELT(grad, i)),
1267                *tmp = Calloc(ncisq, double);
1268    
1269            F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1270            if (j)
1271                error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1272            /* calculate and store grad[i,,] %*% t(R) */
1273            for (j = 0; j < nci; j++) {
1274                F77_CALL(dsymv)("U", &nci, &one, ai, &nci, chol + j, &nci,
1275                                &zero, tmp + j, &nci);
1276            }
1277            /* full symmetric product gives diagonals */
1278            F77_CALL(dtrmm)("R", "U", "T", "N", &nci, &nci, &one, chol, &nci,
1279                            Memcpy(ai, tmp, ncisq), &nci);
1280            /* overwrite lower triangle with gradients for positions in L; */
1281            for (ii = 1; ii < nci; ii++) {
1282                for (j = 0; j < ii; j++) {
1283                    ai[ii + j*nci] = 2. * chol[j*ncip1] * tmp[j + ii*nci];
1284                    ai[j + ii*nci] = 0.;
1285                }
1286            }
1287            Free(chol); Free(tmp);
1288        }
1289    }
1290    
1291    SEXP ssclme_grad(SEXP x, SEXP REMLp, SEXP Uncst)
1292    {
1293        SEXP ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1294    
1295        indicator_gradient(x, asLogical(REMLp), ans);
1296        if (asLogical(Uncst))
1297            unconstrained_gradient(ans, GET_SLOT(x, Matrix_OmegaSym));
1298        UNPROTECT(1);
1299        return ans;
1300    }
1301    
1302    
1303    /**
1304     * Return the gradient of the ML or REML deviance.
1305     *
1306     * @param x pointer to an ssclme object
1307     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1308     * @param Uncp pointer to a logical scalar indicating if the unconstrained parameterization is to be used
1309     *
1310     * @return pointer to a numeric vector of the gradient.
1311     */
1312  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1313  {  {
1314      SEXP      SEXP
1315          Omega = GET_SLOT(x, Matrix_OmegaSym),          Omega = GET_SLOT(x, Matrix_OmegaSym),
1316          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ncsl = GET_SLOT(x, Matrix_ncSym),  
1317          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1318      int      int
1319          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1320          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1321          *nc = INTEGER(ncsl),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1322          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1323          cind, i, n = dims[0],          cind, i, n = dims[0],
1324          nf = length(Omega),          nf = length(Omega),
# Line 1263  Line 1330 
1330          alpha,          alpha,
1331          one = 1.;          one = 1.;
1332      SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));      SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
     int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));  
1333    
1334      ssclme_factor(x);      ssclme_factor(x);
1335      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {
# Line 1298  Line 1364 
1364          F77_CALL(dsyrk)("U", "N", &nci, &ki,          F77_CALL(dsyrk)("U", "N", &nci, &ki,
1365                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1366                          &alpha, tmp, &nci);                          &alpha, tmp, &nci);
1367          alpha = ((double)(REML?(nobs-p):nobs));          alpha = (double)(REML ? (nobs-p) : nobs);
1368          F77_CALL(dsyrk)("U", "N", &nci, &mi,          F77_CALL(dsyrk)("U", "N", &nci, &mi,
1369                          &alpha, b + Gp[i], &nci,                          &alpha, b + Gp[i], &nci,
1370                          &one, tmp, &nci);                          &one, tmp, &nci);
# Line 1316  Line 1382 
1382              int k, odind = cind + nci;              int k, odind = cind + nci;
1383              if (uncst) {              if (uncst) {
1384                  int ione = 1, kk;                  int ione = 1, kk;
1385                  double *rr = Calloc(nci, double);                  double *rr = Calloc(nci, double); /* j'th row of R, the Cholesky factor */
1386                  nlme_symmetrize(tmp, nci);                  nlme_symmetrize(tmp, nci);
1387                  for (j = 0; j < nci; j++, cind++) {                  for (j = 0; j < nci; j++, cind++) {
1388                      for (k = 0; k < nci; k++) rr[k] = chol[j + k*nci];                      for (k = 0; k < j; k++) rr[k] = 0.;
1389                        for (k = j; k < nci; k++) rr[k] = chol[j + k*nci];
1390                      REAL(ans)[cind] = 0.;                      REAL(ans)[cind] = 0.;
1391                      for (k = j; k < nci; k++) {                      for (k = j; k < nci; k++) {
1392                          for (kk = j; kk < nci; kk++) {                          for (kk = j; kk < nci; kk++) {
# Line 1327  Line 1394 
1394                                  tmp[kk * nci + k];                                  tmp[kk * nci + k];
1395                          }                          }
1396                      }                      }
                     for (k = 0; k < nci; k++) rr[k] *= rr[j];  
1397                      for (k = j + 1; k < nci; k++) {                      for (k = j + 1; k < nci; k++) {
1398                          REAL(ans)[odind++] =                          REAL(ans)[odind++] = 2. * rr[j] *
1399                              F77_CALL(ddot)(&nci, rr, &ione, tmp + k, &nci) +                              F77_CALL(ddot)(&nci, rr, &ione, tmp + k*nci, &ione);
                             F77_CALL(ddot)(&nci, rr, &ione,  
                                            tmp + k*nci, &ione);  
1400                      }                      }
1401                  }                  }
1402                  Free(rr);                  Free(rr);
# Line 1352  Line 1416 
1416      return ans;      return ans;
1417  }  }
1418    
1419  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)  /**
1420     * Return the Hessian of the ML or REML deviance.  This is a
1421     * placeholder until I work out the evaluation of the analytic
1422     * Hessian, which probably will involve several helper functions.
1423     *
1424     * @param x pointer to an ssclme object
1425     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1426     * @param Uncp pointer to a logical scalar indicating if the
1427     * unconstrained parameterization is to be used
1428     *
1429     * @return pointer to an approximate Hessian matrix
1430     */
1431    SEXP ssclme_Hessian(SEXP x, SEXP REMLp, SEXP Uncp)
1432    {
1433        int j, ncoef = coef_length(length(GET_SLOT(x, Matrix_OmegaSym)),
1434                                   INTEGER(GET_SLOT(x, Matrix_ncSym))),
1435            unc = asLogical(Uncp);
1436        SEXP ans = PROTECT(allocMatrix(REALSXP, ncoef, ncoef)),
1437            base = PROTECT(unc ? ssclme_coefUnc(x) : ssclme_coef(x)),
1438            current = PROTECT(duplicate(base)),
1439            gradient;
1440    
1441        for (j = 0; j < ncoef; j++) {
1442            double delta = (REAL(base)[j] ? 1.e-7 * REAL(base)[j] : 1.e-7);
1443            int i;
1444    
1445            for (i = 0; i < ncoef; i++) REAL(current)[i] = REAL(base)[i];
1446            REAL(current)[j] += delta/2.;
1447            if (unc) {
1448                ssclme_coefGetsUnc(x, current);
1449            } else {
1450                ssclme_coefGets(x, current);
1451            }
1452            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1453            for (i = 0; i < ncoef; i++) REAL(ans)[j * ncoef + i] = REAL(gradient)[i];
1454            UNPROTECT(1);
1455            REAL(current)[j] -= delta;
1456            if (unc) {
1457                ssclme_coefGetsUnc(x, current);
1458            } else {
1459                ssclme_coefGets(x, current);
1460            }
1461            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1462            for (i = 0; i < ncoef; i++)
1463                REAL(ans)[j * ncoef + i] = (REAL(ans)[j * ncoef + i] - REAL(gradient)[i])/
1464                    delta;
1465            UNPROTECT(1);
1466            /* symmetrize */
1467            for (i = 0; i < j; i++) {
1468                REAL(ans)[j * ncoef + i] = REAL(ans)[i * ncoef + j] =
1469                    (REAL(ans)[j * ncoef + i] + REAL(ans)[i * ncoef + j])/2.;
1470            }
1471        }
1472        UNPROTECT(3);
1473        return ans;
1474    }
1475    
1476    /**
1477     * Calculate and return the fitted values.
1478     *
1479     * @param x pointer to an ssclme object
1480     * @param facs list of grouping factors
1481     * @param mmats list of model matrices
1482     * @param useRf pointer to a logical scalar indicating if the random effects should be used
1483     *
1484     * @return pointer to a numeric array of fitted values
1485     */
1486    SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats, SEXP useRf)
1487  {  {
1488      SEXP val, b;      SEXP val, b;
1489      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
# Line 1374  Line 1505 
1505      } else {      } else {
1506          memset(vv, 0, sizeof(double) * nobs);          memset(vv, 0, sizeof(double) * nobs);
1507      }      }
1508        if (asLogical(useRf)) {
1509      b = PROTECT(ssclme_ranef(x));      b = PROTECT(ssclme_ranef(x));
1510      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1511          int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];          int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];
# Line 1390  Line 1522 
1522              j += nn;              j += nn;
1523          }          }
1524      }      }
1525      UNPROTECT(2);          UNPROTECT(1);
1526        }
1527        UNPROTECT(1);
1528      return val;      return val;
1529  }  }
1530    
# Line 1425  Line 1559 
1559      return Omg;      return Omg;
1560  }  }
1561    
1562    /**
1563     * Copy an ssclme object collapsing the fixed effects slots to the response only.
1564     *
1565     * @param x pointer to an ssclme object
1566     *
1567     * @return a duplicate of x with the fixed effects slots collapsed to the response only
1568     */
1569  SEXP ssclme_collapse(SEXP x)  SEXP ssclme_collapse(SEXP x)
1570  {  {
1571      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),
1572          Omega = GET_SLOT(x, Matrix_OmegaSym),          Omega = GET_SLOT(x, Matrix_OmegaSym),
1573          Dim = GET_SLOT(x, Matrix_DimSym);          Dim = GET_SLOT(x, Matrix_DimSym);
1574      int i, nf = length(Omega), nz = INTEGER(Dim)[1];      int nf = length(Omega), nz = INTEGER(Dim)[1];
     SEXP copy[] = {Matrix_DSym, Matrix_DIsqrtSym, Matrix_DimSym,  
                    Matrix_GpSym, Matrix_LIiSym, Matrix_LIpSym,  
                    Matrix_LIxSym, Matrix_LiSym, Matrix_LpSym,  
                    Matrix_LxSym, Matrix_OmegaSym, Matrix_ParentSym,  
                    Matrix_LIxSym, Matrix_LiSym, Matrix_LpSym,  
                    Matrix_bVarSym, Matrix_devianceSym,  
                    Matrix_devCompSym, Matrix_iSym, Matrix_ncSym,  
                    Matrix_statusSym, Matrix_pSym, Matrix_xSym};  
   
     for (i = 0; i < 23; i++)  
         SET_SLOT(ans, copy[i], duplicate(GET_SLOT(x, copy[i])));  
1575    
1576        slot_dup(ans, x, Matrix_DSym);
1577        slot_dup(ans, x, Matrix_DIsqrtSym);
1578        slot_dup(ans, x, Matrix_DimSym);
1579        slot_dup(ans, x, Matrix_GpSym);
1580        slot_dup(ans, x, Matrix_LiSym);
1581        slot_dup(ans, x, Matrix_LpSym);
1582        slot_dup(ans, x, Matrix_LxSym);
1583        slot_dup(ans, x, Matrix_OmegaSym);
1584        slot_dup(ans, x, Matrix_ParentSym);
1585        slot_dup(ans, x, Matrix_bVarSym);
1586        slot_dup(ans, x, Matrix_devianceSym);
1587        slot_dup(ans, x, Matrix_devCompSym);
1588        slot_dup(ans, x, Matrix_iSym);
1589        slot_dup(ans, x, Matrix_ncSym);
1590        slot_dup(ans, x, Matrix_statusSym);
1591        slot_dup(ans, x, Matrix_pSym);
1592        slot_dup(ans, x, Matrix_xSym);
1593      INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;      INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1594      SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));      SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));
1595      REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;      REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;
# Line 1456  Line 1603 
1603  }  }
1604    
1605    
1606    /**
1607     * Create an lme object from its components.  This is not done by
1608     * new("lme", ...) at the R level because of the possibility of
1609     * causing the copying of very large objects.
1610     *
1611     * @param call Pointer to the original call
1612     * @param facs pointer to the list of grouping factors
1613     * @param x pointer to the model matrices (may be of length zero)
1614     * @param model pointer to the model frame
1615     * @param REML pointer to a logical scalar indicating if REML is used
1616     * @param rep pointer to the converged ssclme object
1617     * @param fitted pointer to the fitted values
1618     * @param residuals pointer to the residuals
1619     *
1620     * @return an lme object
1621     */
1622    SEXP ssclme_to_lme(SEXP call, SEXP facs, SEXP x, SEXP model, SEXP REML,
1623                       SEXP rep, SEXP fitted, SEXP residuals)
1624    {
1625        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));
1626    
1627        SET_SLOT(ans, install("call"), call);
1628        SET_SLOT(ans, install("facs"), facs);
1629        SET_SLOT(ans, Matrix_xSym, x);
1630        SET_SLOT(ans, install("model"), model);
1631        SET_SLOT(ans, install("REML"), REML);
1632        SET_SLOT(ans, install("rep"), rep);
1633        SET_SLOT(ans, install("fitted"), fitted);
1634        SET_SLOT(ans, install("residuals"), residuals);
1635        UNPROTECT(1);
1636        return ans;
1637    }
1638    

Legend:
Removed from v.125  
changed lines
  Added in v.245

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