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 166, Sat May 15 12:59:43 2004 UTC revision 582, Mon Feb 28 18:15:21 2005 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)))
4    
5    /**
6     * Using the sscCrosstab object from the grouping factors, generate
7     * the slots in an ssclme object related to the symmetric sparse
8     * matrix representation of Z'Z.  If the model matrices for the
9     * grouping factors have only one column each then the structure can
10     * be copied, otherwise it must be generated from the sscCrosstab and
11     * the number of columns per grouping factor.
12     *
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     */
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 41  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 73  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++) {
100                          AiOut[pos+1] = AiOut[pos] + 1;                              AiOut[pos] = AiOut[pos-1] + 1;
101                          pos++;                          pos++;
102                      }                      }
103                  }                  }
104                    }
105                  for (jj = 1; jj < nci; jj++) { /* repeat the column adding                  for (jj = 1; jj < nci; jj++) { /* repeat the column adding
106                                                  * another diagonal element */                                                  * another diagonal element */
107                      int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];                      int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];
108                      Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);                      Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);
109                      AiOut[ApOut[mjj + 1] - 1] = mjj; /* maybe mjj-1? */                      AiOut[ApOut[mjj + 1] - 1] = mjj;
110                  }                  }
111              }              }
112          }          }
113          Free(map); Free(ncc);          Free(map); Free(ncc);
114      }      }
     SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));  
     snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));  
     for (i = 0; i <= nf; i++) {  
         snc[i] = nc[i];  
     }  
115  }  }
116    
117  /**  /**
# Line 118  Line 133 
133      Free(sz);      Free(sz);
134  }  }
135    
136    /**
137     * 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     * the number of columns in the model matrices.  There is one more
140     * element in ncv than in facs.  The last element is the number of
141     * columns in the model matrix for the fixed effects plus the
142     * response.  (i.e. p+1)
143     *
144     * @param facs pointer to a list of grouping factors
145     * @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)),
# Line 129  Line 157 
157          *dims = INTEGER(dd);          *dims = INTEGER(dd);
158    
159      if (length(ncv) != (nf + 1))      if (length(ncv) != (nf + 1))
160          error("length of nc (%d) should be length of facs (%d) + 1",          error(_("length of nc (%d) should be length of facs (%d) + 1"),
161                length(ncv), nf);                length(ncv), nf);
162      SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS("ssclme")));      SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS("ssclme")));
163      ssc = VECTOR_ELT(val, 0);      ssc = VECTOR_ELT(val, 0);
# Line 168  Line 196 
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));
199      ldl_symbolic(nzcol, Ap, Ai, Lp, Parent,      R_ldl_symbolic(nzcol, Ap, Ai, Lp, Parent,
                  (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */  
                  (int *) R_alloc(nzcol, sizeof(int)), /* Flag */  
200                   (int *) NULL, (int *) NULL); /* P & Pinv */                   (int *) NULL, (int *) NULL); /* P & Pinv */
201      ssclme_calc_maxod(nzcol, Parent);      ssclme_calc_maxod(nzcol, Parent);
202      Lnz = Lp[nzcol];      Lnz = Lp[nzcol];
# Line 214  Line 240 
240      return val;      return val;
241  }  }
242    
243    /**
244     * Copy information on Z'Z accumulated in the bVar array to Z'Z
245     *
246     * @param ncj number of columns in this block
247     * @param Gpj initial column for this group
248     * @param Gpjp initial column for the next group
249     * @param bVj pointer to the ncj x ncj x mj array to be filled
250     * @param Ap column pointer array for Z'Z
251     * @param Ai row indices for Z'Z
252     * @param Ax elements of Z'Z
253     */
254  static  static
255  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[],
256                const int Ap[], const int Ai[], double Ax[])                const int Ap[], const int Ai[], double Ax[])
# Line 223  Line 260 
260          for (k = 0; k < ncj; k++) {          for (k = 0; k < ncj; k++) {
261              diag = Ap[i + k + 1] - 1;              diag = Ap[i + k + 1] - 1;
262              if (Ai[diag] != i+k)              if (Ai[diag] != i+k)
263                  error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",                  error(_("Expected Ai[%d] to be %d (i.e on diagonal) not %d"),
264                        diag, i+k, Ai[diag]);                        diag, i+k, Ai[diag]);
265              Memcpy(Ax + diag - k, bVj + (i+k-Gpj)*ncj, k + 1);              Memcpy(Ax + diag - k, bVj + (i+k-Gpj)*ncj, k + 1);
266          }          }
267      }      }
268  }  }
269    
270    /**
271     * Copy the dimnames from the list of grouping factors and the model
272     * matrices for the grouping factors into the appropriate parts of the
273     * ssclme object.
274     *
275     * @param x pointer to an ssclme object
276     * @param facs pointer to a list of factors
277     * @param mmats pointer to a list of model matrices
278     *
279     * @return NULL
280     */
281  SEXP  SEXP
282  ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)  ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)
283  {  {
# Line 255  Line 303 
303      return R_NilValue;      return R_NilValue;
304  }  }
305    
306    /**
307     * Update the numerical entries x, ZtX, and XtX in an ssclme object
308     * according to a set of model matrices.
309     *
310     * @param x pointer to an ssclme object
311     * @param facs pointer to a list of grouping factors
312     * @param mmats pointer to a list of model matrices
313     *
314     * @return NULL
315     */
316  SEXP  SEXP
317  ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)  ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)
318  {  {
# Line 284  Line 342 
342          int *dims = INTEGER(getAttrib(VECTOR_ELT(mmats, i), R_DimSymbol)),          int *dims = INTEGER(getAttrib(VECTOR_ELT(mmats, i), R_DimSymbol)),
343              nci = nc[i];              nci = nc[i];
344          if (nobs != dims[0])          if (nobs != dims[0])
345              error("Expected %d rows in the %d'th model matrix. Got %d",              error(_("Expected %d rows in the %d'th model matrix. Got %d"),
346                    nobs, i+1, dims[0]);                    nobs, i+1, dims[0]);
347          if (nci != dims[1])          if (nci != dims[1])
348              error("Expected %d columns in the %d'th model matrix. Got %d",              error(_("Expected %d columns in the %d'th model matrix. Got %d"),
349                    nci, i+1, dims[1]);                    nci, i+1, dims[1]);
350          Z[i] = REAL(VECTOR_ELT(mmats, i));          Z[i] = REAL(VECTOR_ELT(mmats, i));
351      }      }
# Line 312  Line 370 
370              int fpji = fpj[i] - 1, /* factor indices are 1-based */              int fpji = fpj[i] - 1, /* factor indices are 1-based */
371                  dind = Ap[Gp[j] + fpji * ncj + 1] - 1;                  dind = Ap[Gp[j] + fpji * ncj + 1] - 1;
372              if (Ai[dind] != (Gp[j] + fpji * ncj))              if (Ai[dind] != (Gp[j] + fpji * ncj))
373                  error("logic error in ssclme_update_mm");                  error(_("logic error in ssclme_update_mm"));
374              if (Ncj) {          /* use bVar to accumulate */              if (Ncj) {          /* use bVar to accumulate */
375                  F77_CALL(dsyrk)("U", "T", &ncj, &ione, &one, Zj+i,                  F77_CALL(dsyrk)("U", "T", &ncj, &ione, &one, Zj+i,
376                                  &nobs, &one, bVj + fpji*ncj*ncj, &ncj);                                  &nobs, &one, bVj + fpji*ncj*ncj, &ncj);
# Line 331  Line 389 
389                  nck = nc[k],                  nck = nc[k],
390                  scalar = ncj == 1 && nck == 1;                  scalar = ncj == 1 && nck == 1;
391              double              double
392                  *Zk = Z[k], *work;                  *Zk = Z[k], *work = (double *) NULL;
393    
394              if (!scalar) work = Calloc(ncj * nck, double);              if (!scalar) work = Calloc(ncj * nck, double);
395              for (i = 0; i < nobs; i++) {              for (i = 0; i < nobs; i++) {
396                  int ii, ind = -1, fpji = fpj[i] - 1,                  int ii, ind = -1, fpji = fpj[i] - 1,
# Line 344  Line 403 
403                          break;                          break;
404                      }                      }
405                  }                  }
406                  if (ind < 0) error("logic error in ssclme_update_mm");                  if (ind < 0) error(_("logic error in ssclme_update_mm"));
407                  if (scalar) {   /* update scalars directly */                  if (scalar) {   /* update scalars directly */
408                      Ax[ind] += Zj[i] * Zk[i];                      Ax[ind] += Zj[i] * Zk[i];
409                  } else {                  } else {
# Line 355  Line 414 
414                      for (jj = 0; jj < nck; jj++) {                      for (jj = 0; jj < nck; jj++) {
415                          ind = Apk[fpki * nck + jj] + offset;                          ind = Apk[fpki * nck + jj] + offset;
416                          if (Ai[ind] != row)                          if (Ai[ind] != row)
417                              error("logic error in ssclme_update_mm");                              error(_("logic error in ssclme_update_mm"));
418                          for (ii = 0; ii < ncj; ii++) {                          for (ii = 0; ii < ncj; ii++) {
419                              Ax[ind++] += work[jj * ncj + ii];                              Ax[ind++] += work[jj * ncj + ii];
420                          }                          }
# Line 371  Line 430 
430      return R_NilValue;      return R_NilValue;
431  }  }
432    
433    /**
434     * Inflate Z'Z according to Omega and create the factorization LDL'
435     *
436     * @param x pointer to an ssclme object
437     *
438     * @return NULL
439     */
440  SEXP ssclme_inflate_and_factor(SEXP x)  SEXP ssclme_inflate_and_factor(SEXP x)
441  {  {
442      SEXP      SEXP
# Line 402  Line 468 
468              for (k = 0; k < ncj; k++) {              for (k = 0; k < ncj; k++) {
469                  diag = Ap[i + k + 1] - 1;                  diag = Ap[i + k + 1] - 1;
470                  if (Ai[diag] != i+k)                  if (Ai[diag] != i+k)
471                      error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",                      error(_("Expected Ai[%d] to be %d (i.e on diagonal) not %d"),
472                            diag, i+k, Ai[diag]);                            diag, i+k, Ai[diag]);
473                  for (ii = 0; ii <= k; ii++) {                  for (ii = 0; ii <= k; ii++) {
474                      xcp[diag + ii - k] += omgj[k*ncj + ii];                      xcp[diag + ii - k] += omgj[k*ncj + ii];
# Line 410  Line 476 
476              }              }
477          }          }
478      }      }
479      j = ldl_numeric(n, Ap, Ai, xcp,      j = R_ldl_numeric(n, Ap, Ai, xcp,
480                      INTEGER(GET_SLOT(x, Matrix_LpSym)),                      INTEGER(GET_SLOT(x, Matrix_LpSym)),
481                      INTEGER(GET_SLOT(x, Matrix_ParentSym)),                      INTEGER(GET_SLOT(x, Matrix_ParentSym)),
482                      Lnz, INTEGER(GET_SLOT(x, Matrix_LiSym)),                      INTEGER(GET_SLOT(x, Matrix_LiSym)),
483                      REAL(GET_SLOT(x, Matrix_LxSym)),                      REAL(GET_SLOT(x, Matrix_LxSym)),
484                      D, Y, Pattern, Flag,                      D, (int *) NULL, (int *) NULL); /* P & Pinv */
                     (int *) NULL, (int *) NULL); /* P & Pinv */  
485      if (j != n)      if (j != n)
486          error("rank deficiency of ZtZ+W detected at column %d",          error(_("rank deficiency of ZtZ+W detected at column %d"),
487                j + 1);                j + 1);
488      for (j = 0; j < n; j++) DIsqrt[j] = 1./sqrt(D[j]);      for (j = 0; j < n; j++) DIsqrt[j] = 1./sqrt(D[j]);
489      Free(Lnz); Free(Flag); Free(Pattern); Free(Y); Free(xcp);      Free(Lnz); Free(Flag); Free(Pattern); Free(Y); Free(xcp);
490      return R_NilValue;      return R_NilValue;
491  }  }
492    
493    
494    /**
495     * If status[["factored"]] is FALSE, create and factor Z'Z+Omega, then
496     * create RZX and RXX, the deviance components, and the value of the
497     * deviance for both ML and REML.
498     *
499     * @param x pointer to an ssclme object
500     *
501     * @return NULL
502     */
503  SEXP ssclme_factor(SEXP x)  SEXP ssclme_factor(SEXP x)
504  {  {
505      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
# Line 477  Line 552 
552                                          nci * nci),                                          nci * nci),
553                                   &nci, &j);                                   &nci, &j);
554                  if (j)                  if (j)
555                      error("Omega[%d] is not positive definite", i + 1);                      error(_("Omega[%d] is not positive definite"), i + 1);
556                  for (j = 0; j < nci; j++) {                  for (j = 0; j < nci; j++) {
557                      accum += 2 * log(tmp[j * (nci + 1)]);                      accum += 2 * log(tmp[j * (nci + 1)]);
558                  }                  }
# Line 490  Line 565 
565          for (i = 0; i < pp1; i++) {          for (i = 0; i < pp1; i++) {
566              int j;              int j;
567              double *RZXi = RZX + i * n;              double *RZXi = RZX + i * n;
568              ldl_lsolve(n, RZXi, Lp, Li, Lx);              R_ldl_lsolve(n, RZXi, Lp, Li, Lx);
569              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
570          }          }
571                                  /* downdate and factor X'X */                                  /* downdate and factor X'X */
# Line 518  Line 593 
593      return R_NilValue;      return R_NilValue;
594  }  }
595    
596    /**
597     * Return the position of probe in the sorted index vector ind.  It is
598     * known that the position is greater than or equal to start so a linear
599     * search from start is used.
600     *
601     * @param probe value to be matched
602     * @param start index at which to start
603     * @param ind vector of indices
604     *
605     * @return index of the entry matching probe
606     */
607  static  static
608  int ldl_update_ind(int probe, int start, const int ind[])  int ldl_update_ind(int probe, int start, const int ind[])
609  {  {
610      while (ind[start] < probe) start++;      while (ind[start] < probe) start++;
611      if (ind[start] > probe) error("logic error in ldl_inverse");      if (ind[start] > probe) error(_("logic error in ldl_inverse"));
612      return start;      return start;
613  }  }
614    
615  /**  /**
616   * Update the diagonal blocks of the inverse of LDL' (=Z'Z+W)   * Update the diagonal blocks of the inverse of LDL' (=Z'Z+W).  The
617     * lower Cholesky factors of the updated blocks are stored in the bVar
618     * slot.
619   *   *
620   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
621   *   *
# Line 535  Line 623 
623    
624   */   */
625  static  static
626  SEXP ldl_inverse(SEXP x)  void ldl_inverse(SEXP x)
627  {  {
628      SEXP      SEXP
629          Gpsl = GET_SLOT(x, Matrix_GpSym),          Gpsl = GET_SLOT(x, Matrix_GpSym),
# Line 608  Line 696 
696                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,
697                                   &nci, &jj);                                   &nci, &jj);
698                  if (jj)         /* should never happen */                  if (jj)         /* should never happen */
699                      error(                      error(_(
700                          "Rank deficient variance matrix at group %d, level %d, error code %d",                          "Rank deficient variance matrix at group %d, level %d, error code %d"),
701                          i + 1, j + 1, jj);                          i + 1, j + 1, jj);
702              }              }
703              Free(tmp); Free(ind);              Free(tmp); Free(ind);
704          }          }
705      }      }
     return R_NilValue;  
706  }  }
707    
708    /**
709     * If necessary, factor Z'Z+Omega, ZtX, and XtX then, if necessary,
710     * form RZX, RXX, and bVar for the inverse of the Cholesky factor.
711     *
712     * @param x pointer to an ssclme object
713     *
714     * @return NULL (x is updated in place)
715     */
716  SEXP ssclme_invert(SEXP x)  SEXP ssclme_invert(SEXP x)
717  {  {
718      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
719      if (!status[0]) ssclme_factor(x);      if (!status[0]) ssclme_factor(x);
720      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))
721          error("Unable to invert singular factor of downdated X'X");          error(_("Unable to invert singular factor of downdated X'X"));
722      if (!status[1]) {      if (!status[1]) {
723          SEXP          SEXP
724              RZXsl = GET_SLOT(x, Matrix_RZXSym);              RZXsl = GET_SLOT(x, Matrix_RZXSym);
# Line 643  Line 738 
738    
739          F77_CALL(dtrtri)("U", "N", &pp1, RXX, &pp1, &i);          F77_CALL(dtrtri)("U", "N", &pp1, RXX, &pp1, &i);
740          if (i)          if (i)
741              error("DTRTRI returned error code %d", i);              error(_("DTRTRI returned error code %d"), i);
742          F77_CALL(dtrmm)("R", "U", "N", "N", &n, &pp1, &one,          F77_CALL(dtrmm)("R", "U", "N", "N", &n, &pp1, &one,
743                          RXX, &pp1, RZX, &n);                          RXX, &pp1, RZX, &n);
744          for (i = 0; i < pp1; i++) {          for (i = 0; i < pp1; i++) {
745              int j; double *RZXi = RZX + i * n;              int j; double *RZXi = RZX + i * n;
746              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
747              ldl_ltsolve(n, RZXi, Lp, Li, Lx);              R_ldl_ltsolve(n, RZXi, Lp, Li, Lx);
748          }          }
749          ldl_inverse(x);          ldl_inverse(x);
750          status[1] = 1;          status[1] = 1;
# Line 657  Line 752 
752      return R_NilValue;      return R_NilValue;
753  }  }
754    
755    /**
756     * Create and insert initial values for Omega_i.
757     *
758     * @param x pointer to an ssclme object
759     *
760     * @return NULL
761     */
762  SEXP ssclme_initial(SEXP x)  SEXP ssclme_initial(SEXP x)
763  {  {
764      SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),      SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),
# Line 685  Line 787 
787          for (j = Gpi; j < p2; j += nci) {          for (j = Gpi; j < p2; j += nci) {
788              for (k = 0; k < nci; k++) {              for (k = 0; k < nci; k++) {
789                  int jk = j+k, jj = Ap[jk+1] - 1;                  int jk = j+k, jj = Ap[jk+1] - 1;
790                  if (Ai[jj] != jk) error("malformed ZtZ structure");                  if (Ai[jj] != jk) error(_("malformed ZtZ structure"));
791                  mm[k * ncip1] += Ax[jj] * mi;                  mm[k * ncip1] += Ax[jj] * mi;
792              }              }
793          }          }
# Line 726  Line 828 
828   *   *
829   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
830   *   *
831   * @return a vector containing the conditional modes of the random effects   * @return a list containing the conditional modes of the random effects
832   */   */
833  SEXP ssclme_ranef(SEXP x)  SEXP ssclme_ranef(SEXP x)
834  {  {
# Line 781  Line 883 
883                                          nobs + 1 - pp1 : nobs))));                                          nobs + 1 - pp1 : nobs))));
884  }  }
885    
886    /**
887     * Calculate the length of the parameter vector, which is called coef
888     * for historical reasons.
889     *
890     * @param nf number of factors
891     * @param nc number of columns in the model matrices for each factor
892     *
893     * @return total length of the coefficient vector
894     */
895  static  static
896  int coef_length(int nf, const int nc[])  int coef_length(int nf, const int nc[])
897  {  {
# Line 790  Line 901 
901  }  }
902    
903  /**  /**
904   * Extract the upper triangles of the Omega matrices.   * Extract the upper triangles of the Omega matrices.  These aren't
905   * (These aren't "coefficients" but the extractor is   * "coefficients" but the extractor is called coef for historical
906   * called coef for historical reasons.)   * reasons.  Within each group these values are in the order of the
907     * diagonal entries first then the strict upper triangle in row
908     * order.
909   *   *
910   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
911     * @param Unconstr pointer to a logical object indicating if the
912     *        unconstrained parameterization should be used.
913   *   *
914   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of the values in the upper triangles of the
915   * Omega matrices   * Omega matrices
916   */   */
917  SEXP ssclme_coef(SEXP x)  SEXP ssclme_coef(SEXP x, SEXP Unconstr)
918  {  {
919      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
920      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
921          i, nf = length(Omega), vind;          i, nf = length(Omega), unc = asLogical(Unconstr), vind;
922      SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));      SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
923      double *vv = REAL(val);      double *vv = REAL(val);
924    
925      vind = 0;      vind = 0;                   /* index in vv */
926      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
927          int nci = nc[i];          int nci = nc[i], ncip1 = nci + 1;
928          if (nci == 1) {          if (nci == 1) {
929              vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];              vv[vind++] = (unc ?
930                              log(REAL(VECTOR_ELT(Omega, i))[0]) :
931                              REAL(VECTOR_ELT(Omega, i))[0]);
932          } else {          } else {
933              int j, k, odind = vind + nci, ncip1 = nci + 1;              if (unc) {          /* L log(D) L' factor of Omega[i,,] */
934                    int j, k, ncisq = nci * nci;
935                    double *tmp = Memcpy(Calloc(ncisq, double),
936                                         REAL(VECTOR_ELT(Omega, i)), ncisq);
937                    F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
938                    if (j)          /* should never happen */
939                        error(_("DPOTRF returned error code %d on Omega[[%d]]"),
940                              j, i+1);
941                    for (j = 0; j < nci; j++) {
942                        double diagj = tmp[j * ncip1];
943                        vv[vind++] = 2. * log(diagj);
944                        for (k = j + 1; k < nci; k++) {
945                            tmp[j + k * nci] /= diagj;
946                        }
947                    }
948                    for (j = 0; j < nci; j++) {
949                        for (k = j + 1; k < nci; k++) {
950                            vv[vind++] = tmp[j + k * nci];
951                        }
952                    }
953                    Free(tmp);
954                } else {            /* upper triangle of Omega[i,,] */
955                    int j, k, odind = vind + nci;
956              double *omgi = REAL(VECTOR_ELT(Omega, i));              double *omgi = REAL(VECTOR_ELT(Omega, i));
957    
958              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
# Line 825  Line 964 
964              vind = odind;              vind = odind;
965          }          }
966      }      }
967        }
968      UNPROTECT(1);      UNPROTECT(1);
969      return val;      return val;
970  }  }
971    
972  /**  /**
973   * Extract the unconstrained parameters that determine the   * Extract the unconstrained parameters that determine the
974   * Omega matrices. (Called coef for historical reasons.)   * Omega matrices. (Called coef for historical reasons.)  The
975     * unconstrained parameters are derived from the LDL' decomposition of
976     * Omega_i.  The first nc[i] entries in each group are the diagonals
977     * of log(D) followed by the strict lower triangle of L in column
978     * order.
979   *   *
980   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
981   *   *
# Line 857  Line 1001 
1001                                   REAL(VECTOR_ELT(Omega, i)), ncisq);                                   REAL(VECTOR_ELT(Omega, i)), ncisq);
1002              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1003              if (j)              /* should never happen */              if (j)              /* should never happen */
1004                  error("DPOTRF returned error code %d on Omega[[%d]]",                  error(_("DPOTRF returned error code %d on Omega[[%d]]"),
1005                        j, i+1);                        j, i+1);
1006              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1007                  double diagj = tmp[j * ncip1];                  double diagj = tmp[j * ncip1];
# Line 895  Line 1039 
1039      double *cc = REAL(coef);      double *cc = REAL(coef);
1040    
1041      if (length(coef) != coef_length(nf, nc) || !isReal(coef))      if (length(coef) != coef_length(nf, nc) || !isReal(coef))
1042          error("coef must be a numeric vector of length %d",          error(_("coef must be a numeric vector of length %d"),
1043                coef_length(nf, nc));                coef_length(nf, nc));
1044      cind = 0;      cind = 0;
1045      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
# Line 935  Line 1079 
1079   *   *
1080   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1081   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
1082     * @param Unc pointer to a logical object indicating if the
1083     *        unconstrained parameterization should be used.
1084   *   *
1085   * @return R_NilValue   * @return R_NilValue
1086   */   */
1087  SEXP ssclme_coefGets(SEXP x, SEXP coef)  SEXP ssclme_coefGets(SEXP x, SEXP coef, SEXP Unc)
1088  {  {
1089      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1090      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1091            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1092          cind, i, nf = length(Omega),          cind, i, nf = length(Omega),
1093          *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));          unc = asLogical(Unc);
1094      double *cc = REAL(coef);      double *cc = REAL(coef);
1095    
1096      if (length(coef) != coef_length(nf, nc) || !isReal(coef))      if (length(coef) != coef_length(nf, nc) || !isReal(coef))
1097          error("coef must be a numeric vector of length %d",          error(_("coef must be a numeric vector of length %d"),
1098                coef_length(nf, nc));                coef_length(nf, nc));
1099      cind = 0;      cind = 0;
1100      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1101          int nci = nc[i];          int nci = nc[i];
1102          if (nci == 1) {          if (nci == 1) {
1103              REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];              REAL(VECTOR_ELT(Omega, i))[0] = (unc ?
1104                                                 exp(cc[cind++]) :
1105                                                 cc[cind++]);
1106          } else {          } else {
1107              int j, k, odind = cind + nci, ncip1 = nci + 1;              int odind = cind + nci, /* off-diagonal index */
1108              double *omgi = REAL(VECTOR_ELT(Omega, i));                  j, k,
1109                    ncip1 = nci + 1,
1110                    ncisq = nci * nci;
1111                double
1112                    *omgi = REAL(VECTOR_ELT(Omega, i));
1113                if (unc) {
1114                    double
1115                        *tmp = Calloc(ncisq, double),
1116                        diagj, one = 1., zero = 0.;
1117    
1118                    memset(omgi, 0, sizeof(double) * ncisq);
1119                    for (j = 0; j < nci; j++) {
1120                        tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1121                        for (k = j + 1; k < nci; k++) {
1122                            tmp[k*nci + j] = cc[odind++] * diagj;
1123                        }
1124                    }
1125                    F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1126                                    tmp, &nci, &zero, omgi, &nci);
1127                    Free(tmp);
1128                } else {
1129              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1130                  omgi[j * ncip1] = cc[cind++];                  omgi[j * ncip1] = cc[cind++];
1131                  for (k = j + 1; k < nci; k++) {                  for (k = j + 1; k < nci; k++) {
1132                      omgi[k*nci + j] = cc[odind++];                      omgi[k*nci + j] = cc[odind++];
1133                  }                  }
1134              }              }
1135                }
1136              cind = odind;              cind = odind;
1137          }          }
1138      }      }
# Line 971  Line 1140 
1140      return x;      return x;
1141  }  }
1142    
1143  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)  
1144    /**
1145     * Returns the inverse of the updated Omega matrices for an ECME
1146     * iteration.  These matrices are also used in the gradient calculation.
1147     *
1148     * @param x pointer to an ssclme object
1149     * @param REML indicator of REML being used
1150     * @param val pointer to a list of symmetric q_i by q_i matrices
1151     */
1152    static void
1153    common_ECME_gradient(SEXP x, int REML, SEXP val)
1154  {  {
1155      SEXP      SEXP
         Omega = GET_SLOT(x, Matrix_OmegaSym),  
1156          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ncsl = GET_SLOT(x, Matrix_ncSym),  
1157          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1158      int      int
1159          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1160          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1161          *nc = INTEGER(ncsl),          i, j, n = *INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1162          *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);  
1163      double      double
1164          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
1165          *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),          *b = RZX + p * INTEGER(getAttrib(RZXsl, R_DimSymbol))[0],
1166          *b,          one = 1.0, zero = 0.0;
         alpha,  
         one = 1.,  
         zero = 0.;  
   
     p = pp1 - 1;  
     b = RZX + p * n;  
     if (verbose) {  
         SEXP coef = PROTECT(ssclme_coef(x));  
         int lc = length(coef); double *cc = REAL(coef);  
1167    
         Rprintf("  EM iterations\n");  
         Rprintf("%3d %.3f", 0, dev[REML ? 1 : 0]);  
         for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);  
         Rprintf("\n");  
         UNPROTECT(1);  
     }  
     for (iter = 0; iter < nEM; iter++) {  
1168          ssclme_invert(x);          ssclme_invert(x);
1169          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
1170              int ki = Gp[i+1] - Gp[i],              int ki = Gp[i+1] - Gp[i],
1171                  nci = nc[i],                  nci = nc[i],
1172                  mi = ki/nci;                  mi = ki/nci;
1173              double              double
1174                  *vali = REAL(VECTOR_ELT(Omega, i));              *vali = REAL(VECTOR_ELT(val, i)),
1175                alpha = (double)(REML ? (nobs-p) : nobs)/((double) mi);
1176    
             alpha = ((double)(REML?(nobs-p):nobs))/((double)mi);  
1177              F77_CALL(dsyrk)("U", "N", &nci, &mi,              F77_CALL(dsyrk)("U", "N", &nci, &mi,
1178                              &alpha, b + Gp[i], &nci,                              &alpha, b + Gp[i], &nci,
1179                              &zero, vali, &nci);                              &zero, vali, &nci);
# Line 1030  Line 1182 
1182                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1183                              &one, vali, &nci);                              &one, vali, &nci);
1184              if (REML) {              if (REML) {
                 int j;  
1185                  for (j = 0; j < p; j++) {                  for (j = 0; j < p; j++) {
1186                      F77_CALL(dsyrk)("U", "N", &nci, &mi,                      F77_CALL(dsyrk)("U", "N", &nci, &mi,
1187                                  &alpha, RZX + Gp[i] + j*n, &nci,                                  &alpha, RZX + Gp[i] + j*n, &nci,
1188                                  &one, vali, &nci);                                  &one, vali, &nci);
1189                  }                  }
1190              }              }
             F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);  
             if (info)  
                 error("DPOTRF returned error code %d in Omega[[%d]] update",  
                       info, i + 1);  
             F77_CALL(dpotri)("U", &nci, vali, &nci, &info);  
             if (info)  
                 error("DPOTRI returned error code %d in Omega[[%d]] update",  
                       info, i + 1);  
1191          }          }
1192          status[0] = status[1] = 0;  }
1193          if (verbose) {  
1194              SEXP coef = PROTECT(ssclme_coef(x));  /**
1195              int lc = length(coef); double *cc = REAL(coef);   * Print the verbose output in the ECME iterations
1196              Rprintf("%3d %.3f", iter + 1, dev[REML ? 1 : 0]);   *
1197     * @param x pointer to an ssclme object
1198     * @param iter iteration number
1199     * @param REML indicator of whether REML is being used
1200     */
1201    static void EMsteps_verbose_print(SEXP x, int iter, int REML)
1202    {
1203        SEXP coef = PROTECT(ssclme_coef(x, ScalarLogical(0)));
1204        int i, lc = length(coef);
1205        double *cc = REAL(coef), *dev = REAL(GET_SLOT(x, Matrix_devianceSym));
1206    
1207    
1208        ssclme_factor(x);
1209        if (iter == 0) Rprintf("  EM iterations\n");
1210        Rprintf("%3d %.3f", iter, dev[REML ? 1 : 0]);
1211              for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);              for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1212              Rprintf("\n");              Rprintf("\n");
1213              UNPROTECT(1);              UNPROTECT(1);
1214          }          }
1215    
1216    /**
1217     * Perform ECME steps for the REML or ML criterion.
1218     *
1219     * @param x pointer to an ssclme object
1220     * @param nsteps pointer to an integer scalar - the number of ECME steps to perform
1221     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1222     * @param verb pointer to a logical scalar indicating verbose mode
1223     *
1224     * @return NULL
1225     */
1226    SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1227    {
1228        SEXP
1229            Omega = GET_SLOT(x, Matrix_OmegaSym);
1230        int
1231            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1232            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1233            REML = asLogical(REMLp),
1234            i, info, iter,
1235            nEM = asInteger(nsteps),
1236            nf = length(Omega),
1237            verbose = asLogical(verb);
1238    
1239        if (verbose) EMsteps_verbose_print(x, 0, REML);
1240        for (iter = 0; iter < nEM; iter++) {
1241            common_ECME_gradient(x, REML, Omega);
1242            status[0] = status[1] = 0;
1243            for (i = 0; i < nf; i++) {
1244                double *vali = REAL(VECTOR_ELT(Omega, i));
1245    
1246                F77_CALL(dpotrf)("U", nc + i, vali, nc + i, &info);
1247                if (info)
1248                    error(_("DPOTRF returned error code %d in Omega[[%d]] update"),
1249                          info, i + 1);
1250                F77_CALL(dpotri)("U", nc + i, vali, nc + i, &info);
1251                if (info)
1252                    error(_("DPOTRI returned error code %d in Omega[[%d]] update"),
1253                          info, i + 1);
1254            }
1255            if (verbose) EMsteps_verbose_print(x, iter + 1, REML);
1256      }      }
1257      ssclme_factor(x);      ssclme_factor(x);
1258      return R_NilValue;      return R_NilValue;
1259  }  }
1260    
1261    /**
1262     * Evaluate the gradient with respect to the indicators of the
1263     * positions in the Omega matrices.
1264     *
1265     * @param x Pointer to an ssclme object
1266     * @param REML indicator of whether REML is to be used
1267     * @param val Pointer to an object of the same structure as Omega
1268     */
1269    static void indicator_gradient(SEXP x, int REML, SEXP val)
1270    {
1271        SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1272        int
1273            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1274            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1275            i, nf = length(Omega);
1276    
1277        common_ECME_gradient(x, REML, val);
1278        for (i = 0; i < nf; i++) {
1279            int info, nci = nc[i], mi = (Gp[i+1] - Gp[i])/nci;
1280            double
1281                *work = Memcpy((double *) Calloc(nci * nci, double),
1282                               REAL(VECTOR_ELT(Omega, i)), nci * nci),
1283                alpha = (double) -mi, beta = (double) mi;
1284    
1285            F77_CALL(dpotrf)("U", &nci, work, &nci, &info);
1286            if (info)
1287                error(_("DPOTRF gave error code %d on Omega[[%d]]"), info, i + 1);
1288            F77_CALL(dtrtri)("U", "N", &nci, work, &nci, &info);
1289            if (info)
1290                error(_("DPOTRI gave error code %d on Omega[[%d]]"), info, i + 1);
1291            F77_CALL(dsyrk)("U", "N", &nci, &nci, &alpha, work, &nci,
1292                            &beta, REAL(VECTOR_ELT(val, i)), &nci);
1293            Free(work);
1294        }
1295    }
1296    
1297    /**
1298     * Overwrite a gradient with respect to positions in Omega[[i]] by the
1299     * gradient with respect to the unconstrained parameters.
1300     *
1301     * @param grad pointer to a gradient wrt positions.  Contents are overwritten.
1302     * @param Omega pointer to a list of symmetric matrices (upper storage).
1303     */
1304    static void unconstrained_gradient(SEXP grad, SEXP Omega)
1305    {
1306        int i, ii, j, nf = length(Omega);
1307        double one = 1.0, zero = 0.0;
1308    
1309        for (i = 0; i < nf; i++) {
1310            SEXP Omegai = VECTOR_ELT(Omega, i);
1311            int nci = *INTEGER(getAttrib(Omegai, R_DimSymbol)),
1312                ncisq = nci * nci, ncip1 = nci + 1;
1313            double *chol =
1314                Memcpy(Calloc(ncisq, double), REAL(Omegai), ncisq),
1315                *ai = REAL(VECTOR_ELT(grad, i)),
1316                *tmp = Calloc(ncisq, double);
1317    
1318            F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1319            if (j)
1320                error(_("DPOTRF gave error code %d on Omega[[%d]]"), j, i + 1);
1321            /* calculate and store grad[i,,] %*% t(R) */
1322            for (j = 0; j < nci; j++) {
1323                F77_CALL(dsymv)("U", &nci, &one, ai, &nci, chol + j, &nci,
1324                                &zero, tmp + j, &nci);
1325            }
1326            /* full symmetric product gives diagonals */
1327            F77_CALL(dtrmm)("R", "U", "T", "N", &nci, &nci, &one, chol, &nci,
1328                            Memcpy(ai, tmp, ncisq), &nci);
1329            /* overwrite upper triangle with gradients for positions in L' */
1330            for (ii = 1; ii < nci; ii++) {
1331                for (j = 0; j < ii; j++) {
1332                    ai[j + ii*nci] = 2. * chol[j*ncip1] * tmp[j + ii*nci];
1333                    ai[ii + j*nci] = 0.;
1334                }
1335            }
1336            Free(chol); Free(tmp);
1337        }
1338    }
1339    
1340    /**
1341     * Fills cvec with unlist(lapply(mList, function(el) el[upper.tri(el, strict = FALSE)]))
1342     *
1343     * @param mList pointer to a list of REAL matrices
1344     * @param nc number of columns in each matrix
1345     * @param cvec pointer to REAL vector to receive the result
1346     */
1347    static void upperTriList_to_vector(SEXP mList, int *nc, SEXP cvec)
1348    {
1349        int i, ii, j, nf = length(mList), pos = 0;
1350    
1351        pos = 0;                    /* position in vector */
1352        for (i = 0; i < nf; i++) {
1353            SEXP omgi = VECTOR_ELT(mList, i);
1354            int nci = nc[i];
1355            for (j = 0; j < nci; j++) {
1356                for (ii = 0; ii <= j; ii++) {
1357                    REAL(cvec)[pos++] = REAL(omgi)[ii + j * nci];
1358                }
1359            }
1360        }
1361    }
1362    
1363    /**
1364     * Return the gradient of the ML or REML deviance.
1365     *
1366     * @param x pointer to an ssclme object
1367     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1368     * @param Unc pointer to a logical scalar indicating if the
1369     *        unconstrained parameterization is to be used
1370     * @param OneVector pointer to a logical scalar indicating if result
1371     *        should be a single vector.
1372     *
1373     * @return pointer to the gradient as a list of matrices or as a vector.
1374     */
1375    SEXP ssclme_grad(SEXP x, SEXP REMLp, SEXP Unc, SEXP OneVector)
1376    {
1377        SEXP ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1378    
1379        indicator_gradient(x, asLogical(REMLp), ans);
1380        if (asLogical(Unc))
1381            unconstrained_gradient(ans, GET_SLOT(x, Matrix_OmegaSym));
1382        if (asLogical(OneVector)) {
1383            int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));
1384            SEXP val = PROTECT(allocVector(REALSXP, coef_length(length(ans), nc)));
1385    
1386            upperTriList_to_vector(ans, nc, val);
1387            UNPROTECT(2);
1388            return val;
1389        }
1390        UNPROTECT(1);
1391        return ans;
1392    }
1393    
1394  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1395  {  {
1396      SEXP      SEXP
1397          Omega = GET_SLOT(x, Matrix_OmegaSym),          Omega = GET_SLOT(x, Matrix_OmegaSym),
1398          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ncsl = GET_SLOT(x, Matrix_ncSym),  
1399          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1400      int      int
1401          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1402          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1403          *nc = INTEGER(ncsl),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1404          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1405          cind, i, n = dims[0],          cind, i, n = dims[0],
1406          nf = length(Omega),          nf = length(Omega),
# Line 1107  Line 1437 
1437    
1438          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1439          if (j)          if (j)
1440              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);              error(_("DPOTRF gave error code %d on Omega[[%d]]"), j, i + 1);
1441          Memcpy(tmp, chol, ncisq);          Memcpy(tmp, chol, ncisq);
1442          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1443          if (j)          if (j)
1444              error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);              error(_("DPOTRI gave error code %d on Omega[[%d]]"), j, i + 1);
1445          alpha = (double) -mi;          alpha = (double) -mi;
1446          F77_CALL(dsyrk)("U", "N", &nci, &ki,          F77_CALL(dsyrk)("U", "N", &nci, &ki,
1447                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1448                          &alpha, tmp, &nci);                          &alpha, tmp, &nci);
1449          alpha = ((double)(REML?(nobs-p):nobs));          alpha = (double)(REML ? (nobs-p) : nobs);
1450          F77_CALL(dsyrk)("U", "N", &nci, &mi,          F77_CALL(dsyrk)("U", "N", &nci, &mi,
1451                          &alpha, b + Gp[i], &nci,                          &alpha, b + Gp[i], &nci,
1452                          &one, tmp, &nci);                          &one, tmp, &nci);
# Line 1134  Line 1464 
1464              int k, odind = cind + nci;              int k, odind = cind + nci;
1465              if (uncst) {              if (uncst) {
1466                  int ione = 1, kk;                  int ione = 1, kk;
1467                  double *rr = Calloc(nci, double);                  double *rr = Calloc(nci, double); /* j'th row of R, the Cholesky factor */
1468                  nlme_symmetrize(tmp, nci);                  nlme_symmetrize(tmp, nci);
1469                  for (j = 0; j < nci; j++, cind++) {                  for (j = 0; j < nci; j++, cind++) {
1470                      for (k = 0; k < nci; k++) rr[k] = chol[j + k*nci];                      for (k = 0; k < j; k++) rr[k] = 0.;
1471                        for (k = j; k < nci; k++) rr[k] = chol[j + k*nci];
1472                      REAL(ans)[cind] = 0.;                      REAL(ans)[cind] = 0.;
1473                      for (k = j; k < nci; k++) {                      for (k = j; k < nci; k++) {
1474                          for (kk = j; kk < nci; kk++) {                          for (kk = j; kk < nci; kk++) {
# Line 1145  Line 1476 
1476                                  tmp[kk * nci + k];                                  tmp[kk * nci + k];
1477                          }                          }
1478                      }                      }
                     for (k = 0; k < nci; k++) rr[k] *= rr[j];  
1479                      for (k = j + 1; k < nci; k++) {                      for (k = j + 1; k < nci; k++) {
1480                          REAL(ans)[odind++] =                          REAL(ans)[odind++] = 2. * rr[j] *
1481                              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);  
1482                      }                      }
1483                  }                  }
1484                  Free(rr);                  Free(rr);
# Line 1170  Line 1498 
1498      return ans;      return ans;
1499  }  }
1500    
1501    /**
1502     * Return the Hessian of the ML or REML deviance.  This is a
1503     * placeholder until I work out the evaluation of the analytic
1504     * Hessian, which probably will involve several helper functions.
1505     *
1506     * @param x pointer to an ssclme object
1507     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1508     * @param Uncp pointer to a logical scalar indicating if the
1509     * unconstrained parameterization is to be used
1510     *
1511     * @return pointer to an approximate Hessian matrix
1512     */
1513    SEXP ssclme_Hessian(SEXP x, SEXP REMLp, SEXP Uncp)
1514    {
1515        int j, ncoef = coef_length(length(GET_SLOT(x, Matrix_OmegaSym)),
1516                                   INTEGER(GET_SLOT(x, Matrix_ncSym)))
1517            /* , unc = asLogical(Uncp) */
1518            ;
1519        SEXP ans = PROTECT(allocMatrix(REALSXP, ncoef, ncoef)),
1520            base = PROTECT(ssclme_coef(x, Uncp)),
1521            current = PROTECT(duplicate(base)),
1522            gradient;
1523    
1524        for (j = 0; j < ncoef; j++) {
1525            double delta = (REAL(base)[j] ? 1.e-7 * REAL(base)[j] : 1.e-7);
1526            int i;
1527    
1528            for (i = 0; i < ncoef; i++) REAL(current)[i] = REAL(base)[i];
1529            REAL(current)[j] += delta/2.;
1530            ssclme_coefGets(x, current, Uncp);
1531            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1532            for (i = 0; i < ncoef; i++) REAL(ans)[j * ncoef + i] = REAL(gradient)[i];
1533            UNPROTECT(1);
1534            REAL(current)[j] -= delta;
1535            ssclme_coefGets(x, current, Uncp);
1536            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1537            for (i = 0; i < ncoef; i++)
1538                REAL(ans)[j * ncoef + i] =
1539                    (REAL(ans)[j * ncoef + i] - REAL(gradient)[i])/ delta;
1540            UNPROTECT(1);
1541            for (i = 0; i < j; i++) { /* symmetrize */
1542                REAL(ans)[j * ncoef + i] = REAL(ans)[i * ncoef + j] =
1543                    (REAL(ans)[j * ncoef + i] + REAL(ans)[i * ncoef + j])/2.;
1544            }
1545        }
1546        UNPROTECT(3);
1547        return ans;
1548    }
1549    
1550    /**
1551     * Calculate and return the fitted values.
1552     *
1553     * @param x pointer to an ssclme object
1554     * @param facs list of grouping factors
1555     * @param mmats list of model matrices
1556     * @param useRf pointer to a logical scalar indicating if the random effects should be used
1557     *
1558     * @return pointer to a numeric array of fitted values
1559     */
1560  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats, SEXP useRf)  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats, SEXP useRf)
1561  {  {
1562      SEXP val, b;      SEXP val, b;
# Line 1178  Line 1565 
1565      double *vv, one = 1.0, zero = 0.0;      double *vv, one = 1.0, zero = 0.0;
1566    
1567      if (nf < 1)      if (nf < 1)
1568          error("null factor list passed to ssclme_fitted");          error(_("null factor list passed to ssclme_fitted"));
1569      nobs = length(VECTOR_ELT(facs, 0));      nobs = length(VECTOR_ELT(facs, 0));
1570      val = PROTECT(allocVector(REALSXP, nobs));      val = PROTECT(allocVector(REALSXP, nobs));
1571      vv = REAL(val);      vv = REAL(val);
# Line 1201  Line 1588 
1588              for (j = 0; j < nobs; ) {              for (j = 0; j < nobs; ) {
1589                  int nn = 1, lev = ff[j];                  int nn = 1, lev = ff[j];
1590                  /* check for adjacent rows with same factor level */                  /* check for adjacent rows with same factor level */
1591                  while (ff[j + nn] == lev) nn++;                  while ((j + nn) < nobs && ff[j + nn] == lev) nn++;
1592                  F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,                  F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
1593                                  &one, mm + j, &nobs,                                  &one, mm + j, &nobs,
1594                                  bb + (lev - 1) * nci, &nci,                                  bb + (lev - 1) * nci, &nci,
# Line 1234  Line 1621 
1621    
1622          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1623          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
1624              error("DPOTRF returned error code %d on Omega[%d]",              error(_("DPOTRF returned error code %d on Omega[%d]"),
1625                    j, i + 1);                    j, i + 1);
1626          F77_CALL(dpotri)("U", &nci, mm, &nci, &j);          F77_CALL(dpotri)("U", &nci, mm, &nci, &j);
1627          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
1628              error("DTRTRI returned error code %d on Omega[%d]",              error(_("DTRTRI returned error code %d on Omega[%d]"),
1629                    j, i + 1);                    j, i + 1);
1630          nlme_symmetrize(mm, nci);          nlme_symmetrize(mm, nci);
1631      }      }
# Line 1246  Line 1633 
1633      return Omg;      return Omg;
1634  }  }
1635    
1636  #define slot_dup(sym)  SET_SLOT(ans, sym, duplicate(GET_SLOT(x, sym)))  /**
1637     * Copy an ssclme object collapsing the fixed effects slots to the response only.
1638     *
1639     * @param x pointer to an ssclme object
1640     *
1641     * @return a duplicate of x with the fixed effects slots collapsed to the response only
1642     */
1643  SEXP ssclme_collapse(SEXP x)  SEXP ssclme_collapse(SEXP x)
1644  {  {
1645      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),
# Line 1255  Line 1647 
1647          Dim = GET_SLOT(x, Matrix_DimSym);          Dim = GET_SLOT(x, Matrix_DimSym);
1648      int nf = length(Omega), nz = INTEGER(Dim)[1];      int nf = length(Omega), nz = INTEGER(Dim)[1];
1649    
1650      slot_dup(Matrix_DSym);      slot_dup(ans, x, Matrix_DSym);
1651      slot_dup(Matrix_DIsqrtSym);      slot_dup(ans, x, Matrix_DIsqrtSym);
1652      slot_dup(Matrix_DimSym);      slot_dup(ans, x, Matrix_DimSym);
1653      slot_dup(Matrix_GpSym);      slot_dup(ans, x, Matrix_GpSym);
1654      slot_dup(Matrix_LiSym);      slot_dup(ans, x, Matrix_LiSym);
1655      slot_dup(Matrix_LpSym);      slot_dup(ans, x, Matrix_LpSym);
1656      slot_dup(Matrix_LxSym);      slot_dup(ans, x, Matrix_LxSym);
1657      slot_dup(Matrix_OmegaSym);      slot_dup(ans, x, Matrix_OmegaSym);
1658      slot_dup(Matrix_ParentSym);      slot_dup(ans, x, Matrix_ParentSym);
1659      slot_dup(Matrix_bVarSym);      slot_dup(ans, x, Matrix_bVarSym);
1660      slot_dup(Matrix_devianceSym);      slot_dup(ans, x, Matrix_devianceSym);
1661      slot_dup(Matrix_devCompSym);      slot_dup(ans, x, Matrix_devCompSym);
1662      slot_dup(Matrix_iSym);      slot_dup(ans, x, Matrix_iSym);
1663      slot_dup(Matrix_ncSym);      slot_dup(ans, x, Matrix_ncSym);
1664      slot_dup(Matrix_statusSym);      slot_dup(ans, x, Matrix_statusSym);
1665      slot_dup(Matrix_pSym);      slot_dup(ans, x, Matrix_pSym);
1666      slot_dup(Matrix_xSym);      slot_dup(ans, x, Matrix_xSym);
1667      INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;      INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1668      SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));      SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));
1669      REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;      REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;
# Line 1284  Line 1676 
1676      return ans;      return ans;
1677  }  }
1678    
1679    
1680    /**
1681     * Create an lme object from its components.  This is not done by
1682     * new("lme", ...) at the R level because of the possibility of
1683     * causing the copying of very large objects.
1684     *
1685     * @param call Pointer to the original call
1686     * @param facs pointer to the list of grouping factors
1687     * @param x pointer to the model matrices (may be of length zero)
1688     * @param model pointer to the model frame
1689     * @param REML pointer to a logical scalar indicating if REML is used
1690     * @param rep pointer to the converged ssclme object
1691     * @param fitted pointer to the fitted values
1692     * @param residuals pointer to the residuals
1693     * @param terms pointer to a terms object (redundant if model is non-empty)
1694     * @param assign pointer to an integer assign vector
1695     *
1696     * @return an lme object
1697     */
1698  SEXP ssclme_to_lme(SEXP call, SEXP facs, SEXP x, SEXP model, SEXP REML,  SEXP ssclme_to_lme(SEXP call, SEXP facs, SEXP x, SEXP model, SEXP REML,
1699                     SEXP rep, SEXP fitted, SEXP residuals)                     SEXP rep, SEXP fitted, SEXP residuals, SEXP terms,
1700                       SEXP assign)
1701  {  {
1702      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));
1703    
# Line 1297  Line 1709 
1709      SET_SLOT(ans, install("rep"), rep);      SET_SLOT(ans, install("rep"), rep);
1710      SET_SLOT(ans, install("fitted"), fitted);      SET_SLOT(ans, install("fitted"), fitted);
1711      SET_SLOT(ans, install("residuals"), residuals);      SET_SLOT(ans, install("residuals"), residuals);
1712        SET_SLOT(ans, install("terms"), terms);
1713        SET_SLOT(ans, install("assign"), assign);
1714      UNPROTECT(1);      UNPROTECT(1);
1715      return ans;      return ans;
1716  }  }
1717    

Legend:
Removed from v.166  
changed lines
  Added in v.582

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