SCM Repository

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

Diff of /pkg/src/ssclme.c

revision 100, Sat Apr 17 17:03:59 2004 UTC revision 386, Mon Dec 13 15:59:24 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     * be copied, otherwise it must be generated from the sscCrosstab and
11     * the number of columns per grouping factor.
12   *   *
13   * @param ctab pointer to a pairwise cross-tabulation object   * @param nf number of factors
14   *   * @param nc vector of length nf+2 with number of columns in model matrices
15   * @return pointer to an integer R vector.   * @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++) {
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  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 252  Line 185
185      SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));      SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));
186      SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));      SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));
187      SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));      SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));
188          /* Zero the symmetric matrices (for cosmetic reasons only). */                                  /* Zero symmetric matrices (cosmetic) */
189      memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,      memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,
190             sizeof(double) * pp1 * pp1);             sizeof(double) * pp1 * pp1);
191      memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,      memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,
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));
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);
202      Lnz = Lp[nzcol];      Lnz = Lp[nzcol];
203      SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));      SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));
204      SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));      SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));
# Line 304  Line 236
236          memset(REAL(VECTOR_ELT(tmp, i)), 0,          memset(REAL(VECTOR_ELT(tmp, i)), 0,
237                 sizeof(double) * nci * nci * mi);                 sizeof(double) * nci * nci * mi);
238      }      }
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);
}
239      UNPROTECT(2);      UNPROTECT(2);
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 335  Line 267
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 360  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 369  Line 322
322          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
323          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
324          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
325            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
326          i, j, k,          i, j, k,
327          ione = 1,          ione = 1,
328          nf = length(mmats) - 1,          nf = length(mmats) - 1,
# Line 432  Line 386
386          for (k = j+1; k < nf; k++) { /* off-diagonals */          for (k = j+1; k < nf; k++) { /* off-diagonals */
387              int *fpk = INTEGER(VECTOR_ELT(facs, k)),              int *fpk = INTEGER(VECTOR_ELT(facs, k)),
388                  *Apk = Ap + Gp[k],                  *Apk = Ap + Gp[k],
389                  nck = nc[k];                  nck = nc[k],
390                    scalar = ncj == 1 && nck == 1;
391              double              double
392                  *Zk = Z[k];                  *Zk = Z[k], *work = (double *) NULL;
393
394                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,
397                      row = Gp[j] + fpji * ncj,                      row = Gp[j] + fpji * ncj,
398                      fpki = fpk[i] - 1,                      fpki = fpk[i] - 1,
399                      lastind = Apk[fpki + 1];                      lastind = Apk[fpki*nck + 1];
400                  for (ii = Apk[fpki]; ii < lastind; ii++) {                  for (ii = Apk[fpki*nck]; ii < lastind; ii++) {
401                      if (Ai[ii] == row) {                      if (Ai[ii] == row) {
402                          ind = ii;                          ind = ii;
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 (Ncj || nck > 1) {                  if (scalar) {   /* update scalars directly */
408                                  /* FIXME: run a loop to update */                      Ax[ind] += Zj[i] * Zk[i];
409                      error("code not yet written");                  } else {
410                  } else {        /* update scalars directly */                      int jj, offset = ind - Apk[fpki * nck];
411                      Ax[ind] += Zj[fpji] * Zk[fpki];                      F77_CALL(dgemm)("T", "N", &ncj, &nck, &ione, &one,
412                                        Zj + i, &nobs, Zk + i, &nobs,
413                                        &zero, work, &ncj);
414                        for (jj = 0; jj < nck; jj++) {
415                            ind = Apk[fpki * nck + jj] + offset;
416                            if (Ai[ind] != row)
417                                error("logic error in ssclme_update_mm");
418                            for (ii = 0; ii < ncj; ii++) {
419                                Ax[ind++] += work[jj * ncj + ii];
420                            }
421                        }
422                  }                  }
423              }              }
424                if (!scalar) Free(work);
425          }          }
426      }      }
427      Free(Z);      Free(Z);
428      ssclme_transfer_dimnames(x, facs, mmats);      ssclme_transfer_dimnames(x, facs, mmats);
429        status[0] = status[1] = 0;
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 501  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);
# Line 516  Line 490
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 581  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 589  Line 573
573          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
574                          RZX, &n, &one, RXX, &pp1);                          RZX, &n, &one, RXX, &pp1);
575          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
576          if (i)          if (i) {
577              error("DPOTRF returned error code %d", i);              warning("Could not factor downdated X'X, code %d", i);
578                dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;
579            } else {
580                                  /* logdet of RXX */                                  /* logdet of RXX */
581          for (i = 0; i < (pp1 - 1); i++)          for (i = 0; i < (pp1 - 1); i++)
582              dcmp[2] += 2 * log(RXX[i*pp2]);              dcmp[2] += 2 * log(RXX[i*pp2]);
# Line 600  Line 586
586              dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));              dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));
587          deviance[1] = dcmp[0] - dcmp[1] + /* REML */          deviance[1] = dcmp[0] - dcmp[1] + /* REML */
588              dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));              dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));
589            }
590          status[0] = 1;          status[0] = 1;
591          status[1] = 0;          status[1] = 0;
592      }      }
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  {  {
# Line 615  Line 613
613  }  }
614
615  /**  /**
616   * 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
617   * of LDL' (=Z'Z+W)   * 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 624  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),
LIisl = GET_SLOT(x, Matrix_LIiSym),
LIpsl = GET_SLOT(x, Matrix_LIpSym),
630          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
631      int *Gp = INTEGER(Gpsl),      int *Gp = INTEGER(Gpsl),
632          *Li,          *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym)),
633          *LIp = INTEGER(LIpsl), *Lp,          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
634          i,          i,
635          nf = length(Gpsl) - 1,          nf = length(Gpsl) - 1,
636          nzc = length(LIpsl) - 1;          nzc = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
637      double      int maxod = Parent[nzc];
638          *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),      double *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym));
*Lx;
639
640      ssclme_factor(x);      ssclme_factor(x);
641      if (LIp[nzc] == 0) {        /* L and LI are the identity */      if (maxod == 0) {           /* L and L^{-1} are the identity */
642          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
643              Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],              Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],
644                     Gp[i+1] - Gp[i]);                     Gp[i+1] - Gp[i]);
645          }          }
646          return R_NilValue;      } else {
647      }          int *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
648      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;
649
650          double *LIx = REAL(GET_SLOT(x, Matrix_LIxSym)),          double one = 1.0, zero = 0.,
651              one = 1., zero = 0.;              *Lx = REAL(GET_SLOT(x, Matrix_LxSym));

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;

nr = -1;
for (j = G1; j < G2; j += nci) {
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 (cosmetic) */
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)("U", "T", &nci, &rr, &one, tmp, &nr,
&zero, bVi + (j - G1) * nci, &nci);
F77_CALL(dpotrf)("U", &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);
}
}
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);
652
653          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
654              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,
655                  m = maxod + nci,                  m = maxod + 1,
656                  *ind = Calloc(m, int);                  *ind = Calloc(m, int), G1 = Gp[i], G2 = Gp[i+1];
657              double              double
658                  *tmp = Calloc(m * nci, double),                  *tmp = Calloc(m * nci, double),
659                  *mpt = REAL(VECTOR_ELT(bVar, i));                  *bVi = REAL(VECTOR_ELT(bVar, i));
660
661              for (j = Gp[i]; j < Gp[i+1]; j += nci) {                                  /* initialize bVi to zero */
662                  memset((void *) tmp, 0, sizeof(double) * m * nci);              memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
663
664                  kk = 0;         /* ind holds indices of non-zeros */              for (j = G1; j < G2; j += nci) {
665                    kk = 0;         /* ind gets indices of non-zeros */
666                  jj = j;         /* in this block of columns */                  jj = j;         /* in this block of columns */
667                  while (jj >= 0) {                  while (jj >= 0) {
668                      ind[kk++] = jj;                      ind[kk++] = jj;
# Line 760  Line 674
674                  for (k = 0; k < nci; k++) {                  for (k = 0; k < nci; k++) {
675                      double *ccol = tmp + k * nr;                      double *ccol = tmp + k * nr;
676
677                      ccol[k] = 1.;                      for (kk = 0; kk < nr; kk++) ccol[kk] = 0.;
678                      kk = k;                      ccol[k] = 1.; /* initialize from unit diagonal */
679                      for (jj = j + k; jj >= 0; jj = Parent[jj]) {                      for (jj = j + k; jj >= 0; jj = Parent[jj]) {
680                          p2 = Lp[jj+1];                          p2 = Lp[jj+1];
681                          pp = kk;                          pp = pj = ldl_update_ind(jj, 0, ind);
682                          for (p = Lp[jj]; p < p2; p++) {                          for (p = Lp[jj]; p < p2; p++) {
683                              pp = ldl_update_ind(Li[p], pp, ind);                              pp = ldl_update_ind(Li[p], pp, ind);
684                              ccol[pp] -= Lx[p] * ccol[kk];                              ccol[pp] -= Lx[p] * ccol[pj];
685                          }                          }
686                      }                      }
687                  }                  }
688                                  /* scale rows */
689                  for (kk = 0; kk < nr; kk++) {                  for (kk = 0; kk < nr; kk++) { /* scale rows */
690                      for (k = 0; k < nci; k++) {                      for (k = 0; k < nci; k++) {
691                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];
692                      }                      }
693                  }                  }
694                  F77_CALL(dsyrk)("U", "T", &nci, &nr, &one, tmp, &nr,                  F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,
695                                  &zero, mpt + (j - Gp[i])*nci, &nci);                                  &zero, bVi + (j - G1)*nci, &nci);
696                  F77_CALL(dpotrf)("U", &nci, mpt + (j - Gp[i])*nci,                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,
697                                   &nci, &info);                                   &nci, &jj);
698                  if (info)       /* should never happen */                  if (jj)         /* should never happen */
699                      error(                      error(
700                          "Rank deficient variance matrix at group %d, level %d",                          "Rank deficient variance matrix at group %d, level %d, error code %d",
701                          i + 1, j + 1);                          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]))
721            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 820  Line 744
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 828  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 867  Line 798
798
799  /**  /**
800   * Extract the conditional estimates of the fixed effects   * Extract the conditional estimates of the fixed effects
801   *   *
802   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
803   *   *
# Line 895  Line 825
825
826  /**  /**
827   * Extract the conditional modes of the random effects.   * Extract the conditional modes of the random effects.
* FIXME: Change the returned value to be a named list of matrices
*        with dimnames.
828   *   *
829   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
830   *   *
# Line 946  Line 874
874  {  {
875      SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);      SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
876      int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],      int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],
877          nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))          nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))[
878          [length(GET_SLOT(x, Matrix_GpSym))];              length(GET_SLOT(x, Matrix_GpSym))];
879
880      ssclme_invert(x);      ssclme_invert(x);
881      return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *      return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *
# Line 955  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     *
894     */
895  static  static
896  int coef_length(int nf, const int nc[])  int coef_length(int nf, const int nc[])
897  {  {
# Line 964  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 are not in any sense "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   *   *
912   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of the values in the upper triangles of the
913   * Omega matrices   * Omega matrices
914   */   */
915  SEXP ssclme_coef(SEXP x)  SEXP ssclme_coef(SEXP x, SEXP Unconstr)
916  {  {
917      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
918      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
919          i, nf = length(Omega), vind;          i, nf = length(Omega), unc = asLogical(Unconstr), vind;
920      SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));      SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
921      double *vv = REAL(val);      double *vv = REAL(val);
922
923      vind = 0;      vind = 0;                   /* index in vv */
924      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
925          int j, k, nci = nc[i];          int nci = nc[i], ncip1 = nci + 1;
926            if (nci == 1) {
927                vv[vind++] = (unc ?
928                              log(REAL(VECTOR_ELT(Omega, i))[0]) :
929                              REAL(VECTOR_ELT(Omega, i))[0]);
930            } else {
931                if (unc) {          /* L log(D) L' factor of Omega[i,,] */
932                    int j, k, ncisq = nci * nci;
933                    double *tmp = Memcpy(Calloc(ncisq, double),
934                                         REAL(VECTOR_ELT(Omega, i)), ncisq);
935                    F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
936                    if (j)          /* should never happen */
937                        error("DPOTRF returned error code %d on Omega[[%d]]",
938                              j, i+1);
939                    for (j = 0; j < nci; j++) {
940                        double diagj = tmp[j * ncip1];
941                        vv[vind++] = 2. * log(diagj);
942                        for (k = j + 1; k < nci; k++) {
943                            tmp[j + k * nci] /= diagj;
944                        }
945                    }
946                    for (j = 0; j < nci; j++) {
947                        for (k = j + 1; k < nci; k++) {
948                            vv[vind++] = tmp[j + k * nci];
949                        }
950                    }
951                    Free(tmp);
952                } else {            /* upper triangle of Omega[i,,] */
953                    int j, k, odind = vind + nci;
954          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
955
956          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
957              for (k = 0; k <= j; k++) {                      vv[vind++] = omgi[j * ncip1];
958                  vv[vind++] = omgi[j*nci + k];                      for (k = j + 1; k < nci; k++) {
959                            vv[odind++] = omgi[k*nci + j];
960                        }
961                    }
962                    vind = odind;
963              }              }
964          }          }
965      }      }
# Line 996  Line 968
968  }  }
969
970  /**  /**
971   * Extract the upper triangles of the Omega matrices in the unconstrained   * Extract the unconstrained parameters that determine the
972   * parameterization.   * Omega matrices. (Called coef for historical reasons.)  The
973   * (These are not in any sense "coefficients" but the extractor is   * unconstrained parameters are derived from the LDL' decomposition of
974   * called coef for historical reasons.)   * Omega_i.  The first nc[i] entries in each group are the diagonals
975     * of log(D) followed by the strict lower triangle of L in column
976     * order.
977   *   *
978   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
979   *   *
980   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of unconstrained parameters that determine the
981   * Omega matrices   * Omega matrices
982   */   */
983  SEXP ssclme_coefUnc(SEXP x)  SEXP ssclme_coefUnc(SEXP x)
# Line 1025  Line 999
999                                   REAL(VECTOR_ELT(Omega, i)), ncisq);                                   REAL(VECTOR_ELT(Omega, i)), ncisq);
1000              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1001              if (j)              /* should never happen */              if (j)              /* should never happen */
1002                  error("DPOTRF returned error code %d", j);                  error("DPOTRF returned error code %d on Omega[[%d]]",
1003                          j, i+1);
1004              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1005                  double diagj = tmp[j * ncip1];                  double diagj = tmp[j * ncip1];
1006                  vv[vind++] = 2. * log(diagj);                  vv[vind++] = 2. * log(diagj);
# Line 1046  Line 1021
1021  }  }
1022
1023  /**  /**
1024   * Assign the upper triangles of the Omega matrices in the unconstrained   * Assign the Omega matrices from the unconstrained parameterization.
* parameterization.
1025   *   *
1026   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1027   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1078  Line 1052
1052              double              double
1053                  *omgi = REAL(VECTOR_ELT(Omega, i)),                  *omgi = REAL(VECTOR_ELT(Omega, i)),
1054                  *tmp = Calloc(ncisq, double),                  *tmp = Calloc(ncisq, double),
1055                  diagj, one = 1.;                  diagj, one = 1., zero = 0.;
1056              /* FIXEME: Change this to use a factor and dsyrk */
/* LD in omgi and L' in tmp */
1057              memset(omgi, 0, sizeof(double) * ncisq);              memset(omgi, 0, sizeof(double) * ncisq);
1058              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1059                  omgi[j * ncip1] = diagj = exp(cc[cind++]);                  tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1060                  for (k = j + 1; k < nci; k++) {                  for (k = j + 1; k < nci; k++) {
1061                      omgi[j*nci + k] = diagj * (tmp[k*nci + j] = cc[odind++]);                      tmp[k*nci + j] = cc[odind++] * diagj;
1062                  }                  }
1063              }              }
1064              F77_CALL(dtrmm)("R", "U", "N", "U", &nci, &nci, &one,              F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1065                              tmp, &nci, omgi, &nci);                              tmp, &nci, &zero, omgi, &nci);
1066              Free(tmp);              Free(tmp);
1067              cind = odind;              cind = odind;
1068          }          }
# Line 1100  Line 1073
1073
1074  /**  /**
1075   * Assign the upper triangles of the Omega matrices.   * Assign the upper triangles of the Omega matrices.
1076   * (These are not in any sense "coefficients" but are   * (Called coef for historical reasons.)
* called coef for historical reasons.)
1077   *   *
1078   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1079   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
1080   *   *
1081   * @return R_NilValue   * @return R_NilValue
1082   */   */
1083  SEXP ssclme_coefGets(SEXP x, SEXP coef)  SEXP ssclme_coefGets(SEXP x, SEXP coef, SEXP Unc)
1084  {  {
1085      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1086      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1087            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1088          cind, i, nf = length(Omega),          cind, i, nf = length(Omega),
1089          *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));          unc = asLogical(Unc);
1090      double *cc = REAL(coef);      double *cc = REAL(coef);
1091
1092      if (length(coef) != coef_length(nf, nc) || !isReal(coef))      if (length(coef) != coef_length(nf, nc) || !isReal(coef))
# Line 1121  Line 1094
1094                coef_length(nf, nc));                coef_length(nf, nc));
1095      cind = 0;      cind = 0;
1096      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1097          int j, k, nci = nc[i];          int nci = nc[i];
1098          double *omgi = REAL(VECTOR_ELT(Omega, i));          if (nci == 1) {
1099                REAL(VECTOR_ELT(Omega, i))[0] = (unc ?
1100                                                 exp(cc[cind++]) :
1101                                                 cc[cind++]);
1102            } else {
1103                int odind = cind + nci, /* off-diagonal index */
1104                    j, k,
1105                    ncip1 = nci + 1,
1106                    ncisq = nci * nci;
1107                double
1108                    *omgi = REAL(VECTOR_ELT(Omega, i));
1109                if (unc) {
1110                    double
1111                        *tmp = Calloc(ncisq, double),
1112                        diagj, one = 1., zero = 0.;
1113
1114                    memset(omgi, 0, sizeof(double) * ncisq);
1115          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
1116              for (k = 0; k <= j; k++) {                      tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1117                  omgi[j*nci + k] = cc[cind++];                      for (k = j + 1; k < nci; k++) {
1118                            tmp[k*nci + j] = cc[odind++] * diagj;
1119                        }
1120                    }
1121                    F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1122                                    tmp, &nci, &zero, omgi, &nci);
1123                    Free(tmp);
1124                } else {
1125                    for (j = 0; j < nci; j++) {
1126                        omgi[j * ncip1] = cc[cind++];
1127                        for (k = j + 1; k < nci; k++) {
1128                            omgi[k*nci + j] = cc[odind++];
1129                        }
1130              }              }
1131          }          }
1132                cind = odind;
1133            }
1134      }      }
1135      status[0] = status[1] = 0;      status[0] = status[1] = 0;
1136      return x;      return x;
1137  }  }
1138
1139  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1140    /**
1141     * Returns the inverse of the updated Omega matrices for an ECME
1142     * iteration.  These matrices are also used in the gradient calculation.
1143     *
1144     * @param x pointer to an ssclme object
1145     * @param REML indicator of REML being used
1146     * @param val pointer to a list of symmetric q_i by q_i matrices
1147     */
1148    static void
1149    common_ECME_gradient(SEXP x, int REML, SEXP val)
1150  {  {
1151      SEXP      SEXP
Omega = GET_SLOT(x, Matrix_OmegaSym),
1152          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
ncsl = GET_SLOT(x, Matrix_ncSym),
1153          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1154      int      int
1155          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1156          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1157          *nc = INTEGER(ncsl),          i, j, n = *INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1158          *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);
1159      double      double
1160          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
1161          *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),          *b = RZX + p * INTEGER(getAttrib(RZXsl, R_DimSymbol))[0],
1162          *b,          one = 1.0, zero = 0.0;
alpha,
one = 1.,
zero = 0.;
1163
p = pp1 - 1;
b = RZX + p * n;
if (verbose) Rprintf("  EM iterations\n");
for (iter = 0; iter <= nEM; iter++) {
1164          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);
}
1165          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
1166              int ki = Gp[i+1] - Gp[i],              int ki = Gp[i+1] - Gp[i],
1167                  nci = nc[i],                  nci = nc[i],
1168                  mi = ki/nci;                  mi = ki/nci;
1169              double              double
1170                  *vali = REAL(VECTOR_ELT(Omega, i));              *vali = REAL(VECTOR_ELT(val, i)),
1171                alpha = (double)(REML ? (nobs-p) : nobs)/((double) mi);
1172
alpha = ((double)(REML?(nobs-p):nobs))/((double)mi);
1173              F77_CALL(dsyrk)("U", "N", &nci, &mi,              F77_CALL(dsyrk)("U", "N", &nci, &mi,
1174                              &alpha, b + Gp[i], &nci,                              &alpha, b + Gp[i], &nci,
1175                              &zero, vali, &nci);                              &zero, vali, &nci);
# Line 1191  Line 1178
1178                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1179                              &one, vali, &nci);                              &one, vali, &nci);
1180              if (REML) {              if (REML) {
int j;
1181                  for (j = 0; j < p; j++) {                  for (j = 0; j < p; j++) {
1182                      F77_CALL(dsyrk)("U", "N", &nci, &mi,                      F77_CALL(dsyrk)("U", "N", &nci, &mi,
1183                                  &alpha, RZX + Gp[i] + j*n, &nci,                                  &alpha, RZX + Gp[i] + j*n, &nci,
1184                                  &one, vali, &nci);                                  &one, vali, &nci);
1185                  }                  }
1186              }              }
F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
if (info) error("DPOTRF returned error code %d", info);
F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
if (info) error("DPOTRI returned error code %d", info);
1187          }          }
1188    }
1189
1190    /**
1191     * Print the verbose output in the ECME iterations
1192     *
1193     * @param x pointer to an ssclme object
1194     * @param iter iteration number
1195     * @param REML indicator of whether REML is being used
1196     */
1197    static void EMsteps_verbose_print(SEXP x, int iter, int REML)
1198    {
1199        SEXP coef = PROTECT(ssclme_coef(x, ScalarLogical(0)));
1200        int i, lc = length(coef);
1201        double *cc = REAL(coef), *dev = REAL(GET_SLOT(x, Matrix_devianceSym));
1202
1203
1204        ssclme_factor(x);
1205        if (iter == 0) Rprintf("  EM iterations\n");
1206        Rprintf("%3d %.3f", iter, dev[REML ? 1 : 0]);
1207        for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1208        Rprintf("\n");
1209        UNPROTECT(1);
1210    }
1211
1212    /**
1213     * Perform ECME steps for the REML or ML criterion.
1214     *
1215     * @param x pointer to an ssclme object
1216     * @param nsteps pointer to an integer scalar - the number of ECME steps to perform
1217     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1218     * @param verb pointer to a logical scalar indicating verbose mode
1219     *
1220     * @return NULL
1221     */
1222    SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1223    {
1224        SEXP
1225            Omega = GET_SLOT(x, Matrix_OmegaSym);
1226        int
1227            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1228            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1229            REML = asLogical(REMLp),
1230            i, info, iter,
1231            nEM = asInteger(nsteps),
1232            nf = length(Omega),
1233            verbose = asLogical(verb);
1234
1235        if (verbose) EMsteps_verbose_print(x, 0, REML);
1236        for (iter = 0; iter < nEM; iter++) {
1238          status[0] = status[1] = 0;          status[0] = status[1] = 0;
1239            for (i = 0; i < nf; i++) {
1240                double *vali = REAL(VECTOR_ELT(Omega, i));
1241
1242                F77_CALL(dpotrf)("U", nc + i, vali, nc + i, &info);
1243                if (info)
1244                    error("DPOTRF returned error code %d in Omega[[%d]] update",
1245                          info, i + 1);
1246                F77_CALL(dpotri)("U", nc + i, vali, nc + i, &info);
1247                if (info)
1248                    error("DPOTRI returned error code %d in Omega[[%d]] update",
1249                          info, i + 1);
1250            }
1251            if (verbose) EMsteps_verbose_print(x, iter + 1, REML);
1252      }      }
1253      ssclme_factor(x);      ssclme_factor(x);
1254      return R_NilValue;      return R_NilValue;
1255  }  }
1256
1257  SEXP ssclme_gradient(SEXP x, SEXP REMLp)  /**
1258     * Evaluate the gradient with respect to the indicators of the
1259     * positions in the Omega matrices.
1260     *
1261     * @param x Pointer to an ssclme object
1262     * @param REML indicator of whether REML is to be used
1263     * @param val Pointer to an object of the same structure as Omega
1264     */
1265    static void indicator_gradient(SEXP x, int REML, SEXP val)
1266    {
1267        SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1268        int
1269            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1270            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1271            i, nf = length(Omega);
1272
1274        for (i = 0; i < nf; i++) {
1275            int info, nci = nc[i], mi = (Gp[i+1] - Gp[i])/nci;
1276            double
1277                *work = Memcpy((double *) Calloc(nci * nci, double),
1278                               REAL(VECTOR_ELT(Omega, i)), nci * nci),
1279                alpha = (double) -mi, beta = (double) mi;
1280
1281            F77_CALL(dpotrf)("U", &nci, work, &nci, &info);
1282            if (info)
1283                error("DPOTRF gave error code %d on Omega[[%d]]", info, i + 1);
1284            F77_CALL(dtrtri)("U", "N", &nci, work, &nci, &info);
1285            if (info)
1286                error("DPOTRI gave error code %d on Omega[[%d]]", info, i + 1);
1287            F77_CALL(dsyrk)("U", "N", &nci, &nci, &alpha, work, &nci,
1288                            &beta, REAL(VECTOR_ELT(val, i)), &nci);
1289            Free(work);
1290        }
1291    }
1292
1293    /**
1294     * Overwrite a gradient with respect to positions in Omega[[i]] by the
1295     * gradient with respect to the unconstrained parameters.
1296     *
1297     * @param grad pointer to a gradient wrt positions.  Contents are overwritten.
1298     * @param Omega pointer to a list of symmetric matrices (upper storage).
1299     */
1301    {
1302        int i, ii, j, nf = length(Omega);
1303        double one = 1.0, zero = 0.0;
1304
1305        for (i = 0; i < nf; i++) {
1306            SEXP Omegai = VECTOR_ELT(Omega, i);
1307            int nci = *INTEGER(getAttrib(Omegai, R_DimSymbol)),
1308                ncisq = nci * nci, ncip1 = nci + 1;
1309            double *chol =
1310                Memcpy(Calloc(ncisq, double), REAL(Omegai), ncisq),
1312                *tmp = Calloc(ncisq, double);
1313
1314            F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1315            if (j)
1316                error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1317            /* calculate and store grad[i,,] %*% t(R) */
1318            for (j = 0; j < nci; j++) {
1319                F77_CALL(dsymv)("U", &nci, &one, ai, &nci, chol + j, &nci,
1320                                &zero, tmp + j, &nci);
1321            }
1322            /* full symmetric product gives diagonals */
1323            F77_CALL(dtrmm)("R", "U", "T", "N", &nci, &nci, &one, chol, &nci,
1324                            Memcpy(ai, tmp, ncisq), &nci);
1325            /* overwrite upper triangle with gradients for positions in L' */
1326            for (ii = 1; ii < nci; ii++) {
1327                for (j = 0; j < ii; j++) {
1328                    ai[j + ii*nci] = 2. * chol[j*ncip1] * tmp[j + ii*nci];
1329                    ai[ii + j*nci] = 0.;
1330                }
1331            }
1332            Free(chol); Free(tmp);
1333        }
1334    }
1335
1336    /**
1337     * Fills cvec with unlist(lapply(mList, function(el) el[upper.tri(el, strict = FALSE)]))
1338     *
1339     * @param mList pointer to a list of REAL matrices
1340     * @param nc number of columns in each matrix
1341     * @param cvec pointer to REAL vector to receive the result
1342     */
1343    static void upperTriList_to_vector(SEXP mList, int *nc, SEXP cvec)
1344    {
1345        int i, ii, j, nf = length(mList), pos = 0;
1346
1347        pos = 0;                    /* position in vector */
1348        for (i = 0; i < nf; i++) {
1349            SEXP omgi = VECTOR_ELT(mList, i);
1350            int nci = nc[i];
1351            for (j = 0; j < nci; j++) {
1352                for (ii = 0; ii <= j; ii++) {
1353                    REAL(cvec)[pos++] = REAL(omgi)[ii + j * nci];
1354                }
1355            }
1356        }
1357    }
1358
1359    /**
1360     * Return the gradient of the ML or REML deviance.
1361     *
1362     * @param x pointer to an ssclme object
1363     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1364     * @param Uncp pointer to a logical scalar indicating if the unconstrained parameterization is to be used
1365     * @param Uncp pointer to a logical scalar indicating if result should be a single vector
1366     *
1367     * @return pointer to the gradient as a list of matrices or as a vector.
1368     */
1369    SEXP ssclme_grad(SEXP x, SEXP REMLp, SEXP Unc, SEXP OneVector)
1370    {
1371        SEXP ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1372
1374        if (asLogical(Unc))
1376        if (asLogical(OneVector)) {
1377            int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));
1378            SEXP val = PROTECT(allocVector(REALSXP, coef_length(length(ans), nc)));
1379
1380            upperTriList_to_vector(ans, nc, val);
1381            UNPROTECT(2);
1382            return val;
1383        }
1384        UNPROTECT(1);
1385        return ans;
1386    }
1387
1388    SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1389  {  {
1390      SEXP      SEXP
1391            Omega = GET_SLOT(x, Matrix_OmegaSym),
1392          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),
ncsl = GET_SLOT(x, Matrix_ncSym),
1393          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1394      int      int
1395          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1396          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1397          *nc = INTEGER(ncsl),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1398          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1399          i, info,          cind, i, n = dims[0],
1400          n = dims[0],          nf = length(Omega),
1401          nf = length(ncsl) - 2,          nobs, p, pp1 = dims[1],
1402          nobs = nc[nf + 1],          uncst = asLogical(Uncp);
p,
pp1 = dims[1];
1403      double      double
1404          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
1405          *b,          *b,
1406          alpha,          alpha,
1407          one = 1.;          one = 1.;
1408        SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1409
1410        ssclme_factor(x);
1411        if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {
1412            int ncoef = coef_length(nf, nc);
1413            for (i = 0; i < ncoef; i++) REAL(ans)[i] = NA_REAL;
1414            UNPROTECT(1);
1415            return ans;
1416        }
1417        nobs = nc[nf + 1];
1418      p = pp1 - 1;      p = pp1 - 1;
1419      b = RZX + p * n;      b = RZX + p * n;
1420      ssclme_invert(x);      ssclme_invert(x);
1421        cind = 0;
1422      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1423          int ki = Gp[i+1] - Gp[i],          int j, ki = Gp[i+1] - Gp[i],
1424              nci = nc[i],              nci = nc[i], ncip1 = nci + 1, ncisq = nci * nci,
1425              mi = ki/nci;              mi = ki/nci;
1426          double          double
1427              *vali = REAL(VECTOR_ELT(ans, i));              *chol = Memcpy(Calloc(ncisq, double),
1428                               REAL(VECTOR_ELT(Omega, i)), ncisq),
1429                *tmp = Calloc(ncisq, double);
1430
1431          F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1432          if (info)          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1433              error("DPOTRF returned error code %d for component %d of Omega",          if (j)
1434                    info, i + 1);              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1435          F77_CALL(dpotri)("U", &nci, vali, &nci, &info);          Memcpy(tmp, chol, ncisq);
1436          if (info)          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1437              error("DPOTRI returned error code %d for component %d of Omega",          if (j)
1438                    info, i + 1);              error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
1439          alpha = (double) -mi;          alpha = (double) -mi;
1440          F77_CALL(dsyrk)("U", "N", &nci, &ki,          F77_CALL(dsyrk)("U", "N", &nci, &ki,
1441                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1442                          &alpha, vali, &nci);                          &alpha, tmp, &nci);
1443          alpha = ((double)(REML?(nobs-p):nobs));          alpha = (double)(REML ? (nobs-p) : nobs);
1444          F77_CALL(dsyrk)("U", "N", &nci, &mi,          F77_CALL(dsyrk)("U", "N", &nci, &mi,
1445                          &alpha, b + Gp[i], &nci,                          &alpha, b + Gp[i], &nci,
1446                          &one, vali, &nci);                          &one, tmp, &nci);
1447          if (REML) {          if (REML) {
int j;
1448              for (j = 0; j < p; j++) {              for (j = 0; j < p; j++) {
1449                  F77_CALL(dsyrk)("U", "N", &nci, &mi,                  F77_CALL(dsyrk)("U", "N", &nci, &mi,
1450                                  &one, RZX + Gp[i] + j*n, &nci,                                  &one, RZX + Gp[i] + j*n, &nci,
1451                                  &one, vali, &nci);                                  &one, tmp, &nci);
1452                }
1453            }
1454            if (nci == 1) {
1455                REAL(ans)[cind++] = *tmp *
1456                    (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1457            } else {
1458                int k, odind = cind + nci;
1459                if (uncst) {
1460                    int ione = 1, kk;
1461                    double *rr = Calloc(nci, double); /* j'th row of R, the Cholesky factor */
1462                    nlme_symmetrize(tmp, nci);
1463                    for (j = 0; j < nci; j++, cind++) {
1464                        for (k = 0; k < j; k++) rr[k] = 0.;
1465                        for (k = j; k < nci; k++) rr[k] = chol[j + k*nci];
1466                        REAL(ans)[cind] = 0.;
1467                        for (k = j; k < nci; k++) {
1468                            for (kk = j; kk < nci; kk++) {
1469                                REAL(ans)[cind] += rr[k] * rr[kk] *
1470                                    tmp[kk * nci + k];
1471                            }
1472                        }
1473                        for (k = j + 1; k < nci; k++) {
1474                            REAL(ans)[odind++] = 2. * rr[j] *
1475                                F77_CALL(ddot)(&nci, rr, &ione, tmp + k*nci, &ione);
1476              }              }
1477          }          }
1478                    Free(rr);
1479                } else {
1480                    for (j = 0; j < nci; j++) {
1481                        REAL(ans)[cind++] = tmp[j * ncip1];
1482                        for (k = j + 1; k < nci; k++) {
1483                            REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1484                        }
1485                    }
1486                }
1487                cind = odind;
1488            }
1489            Free(tmp); Free(chol);
1490      }      }
1491      UNPROTECT(1);      UNPROTECT(1);
1492      return ans;      return ans;
1493  }  }
1494
1495  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)  /**
1496     * Return the Hessian of the ML or REML deviance.  This is a
1497     * placeholder until I work out the evaluation of the analytic
1498     * Hessian, which probably will involve several helper functions.
1499     *
1500     * @param x pointer to an ssclme object
1501     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1502     * @param Uncp pointer to a logical scalar indicating if the
1503     * unconstrained parameterization is to be used
1504     *
1505     * @return pointer to an approximate Hessian matrix
1506     */
1507    SEXP ssclme_Hessian(SEXP x, SEXP REMLp, SEXP Uncp)
1508    {
1509        int j, ncoef = coef_length(length(GET_SLOT(x, Matrix_OmegaSym)),
1510                                   INTEGER(GET_SLOT(x, Matrix_ncSym)))
1511            /* , unc = asLogical(Uncp) */
1512            ;
1513        SEXP ans = PROTECT(allocMatrix(REALSXP, ncoef, ncoef)),
1514            base = PROTECT(ssclme_coef(x, Uncp)),
1515            current = PROTECT(duplicate(base)),
1517
1518        for (j = 0; j < ncoef; j++) {
1519            double delta = (REAL(base)[j] ? 1.e-7 * REAL(base)[j] : 1.e-7);
1520            int i;
1521
1522            for (i = 0; i < ncoef; i++) REAL(current)[i] = REAL(base)[i];
1523            REAL(current)[j] += delta/2.;
1524            ssclme_coefGets(x, current, Uncp);
1526            for (i = 0; i < ncoef; i++) REAL(ans)[j * ncoef + i] = REAL(gradient)[i];
1527            UNPROTECT(1);
1528            REAL(current)[j] -= delta;
1529            ssclme_coefGets(x, current, Uncp);
1531            for (i = 0; i < ncoef; i++)
1532                REAL(ans)[j * ncoef + i] =
1533                    (REAL(ans)[j * ncoef + i] - REAL(gradient)[i])/ delta;
1534            UNPROTECT(1);
1535            for (i = 0; i < j; i++) { /* symmetrize */
1536                REAL(ans)[j * ncoef + i] = REAL(ans)[i * ncoef + j] =
1537                    (REAL(ans)[j * ncoef + i] + REAL(ans)[i * ncoef + j])/2.;
1538            }
1539        }
1540        UNPROTECT(3);
1541        return ans;
1542    }
1543
1544    /**
1545     * Calculate and return the fitted values.
1546     *
1547     * @param x pointer to an ssclme object
1548     * @param facs list of grouping factors
1549     * @param mmats list of model matrices
1550     * @param useRf pointer to a logical scalar indicating if the random effects should be used
1551     *
1552     * @return pointer to a numeric array of fitted values
1553     */
1554    SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats, SEXP useRf)
1555  {  {
1556      SEXP val, b;      SEXP val, b;
1557      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
# Line 1294  Line 1573
1573      } else {      } else {
1574          memset(vv, 0, sizeof(double) * nobs);          memset(vv, 0, sizeof(double) * nobs);
1575      }      }
1576        if (asLogical(useRf)) {
1577      b = PROTECT(ssclme_ranef(x));      b = PROTECT(ssclme_ranef(x));
1578      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1579          int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];          int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];
# Line 1302  Line 1582
1582          for (j = 0; j < nobs; ) {          for (j = 0; j < nobs; ) {
1583              int nn = 1, lev = ff[j];              int nn = 1, lev = ff[j];
1584              /* check for adjacent rows with same factor level */              /* check for adjacent rows with same factor level */
1585              while (ff[j + nn] == lev) nn++;                  while ((j + nn) < nobs && ff[j + nn] == lev) nn++;
1586              F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,              F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
1587                              &one, mm + j, &nobs,                              &one, mm + j, &nobs,
1588                              bb + (lev - 1) * nci, &nci,                              bb + (lev - 1) * nci, &nci,
# Line 1310  Line 1590
1590              j += nn;              j += nn;
1591          }          }
1592      }      }
1593      UNPROTECT(2);          UNPROTECT(1);
1594        }
1595        UNPROTECT(1);
1596      return val;      return val;
1597  }  }
1598
1599  SEXP ssclme_variances(SEXP x, SEXP REML)  /**
1600     * Return the unscaled variances
1601     *
1602     * @param x pointer to an ssclme object
1603     *
1604     * @return a list similar to the Omega list with the unscaled variances
1605     */
1606    SEXP ssclme_variances(SEXP x)
1607  {  {
1608      SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),      SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
val = PROTECT(allocVector(VECSXP, 2));
1609      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1610          i, nf = length(Omg);          i, nf = length(Omg);
double sigmasq;
1611

SET_VECTOR_ELT(val, 0, Omg);
SET_VECTOR_ELT(val, 1, ssclme_sigma(x, REML));
sigmasq = REAL(VECTOR_ELT(val, 1))[0];
sigmasq = (sigmasq) * (sigmasq);
1612      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1613          double *mm = REAL(VECTOR_ELT(Omg, i));          double *mm = REAL(VECTOR_ELT(Omg, i));
1614          int j, k, nci = nc[i], ncip1 = nci+1;          int j, nci = nc[i];
1615
1616          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1617          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
# Line 1339  Line 1621
1621          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
1622              error("DTRTRI returned error code %d on Omega[%d]",              error("DTRTRI returned error code %d on Omega[%d]",
1623                    j, i + 1);                    j, i + 1);
1624          for (j = 0; j < nci; j++) {          nlme_symmetrize(mm, nci);
mm[j * ncip1] *= sigmasq;
for (k = 0; k < j; k++) {
mm[j + k * nci] = (mm[k + j * nci] *= sigmasq);
1625              }              }
1626        UNPROTECT(1);
1627        return Omg;
1628          }          }
1629
1630    /**
1631     * Copy an ssclme object collapsing the fixed effects slots to the response only.
1632     *
1633     * @param x pointer to an ssclme object
1634     *
1635     * @return a duplicate of x with the fixed effects slots collapsed to the response only
1636     */
1637    SEXP ssclme_collapse(SEXP x)
1638    {
1639        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),
1640            Omega = GET_SLOT(x, Matrix_OmegaSym),
1641            Dim = GET_SLOT(x, Matrix_DimSym);
1642        int nf = length(Omega), nz = INTEGER(Dim)[1];
1643
1644        slot_dup(ans, x, Matrix_DSym);
1645        slot_dup(ans, x, Matrix_DIsqrtSym);
1646        slot_dup(ans, x, Matrix_DimSym);
1647        slot_dup(ans, x, Matrix_GpSym);
1648        slot_dup(ans, x, Matrix_LiSym);
1649        slot_dup(ans, x, Matrix_LpSym);
1650        slot_dup(ans, x, Matrix_LxSym);
1651        slot_dup(ans, x, Matrix_OmegaSym);
1652        slot_dup(ans, x, Matrix_ParentSym);
1653        slot_dup(ans, x, Matrix_bVarSym);
1654        slot_dup(ans, x, Matrix_devianceSym);
1655        slot_dup(ans, x, Matrix_devCompSym);
1656        slot_dup(ans, x, Matrix_iSym);
1657        slot_dup(ans, x, Matrix_ncSym);
1658        slot_dup(ans, x, Matrix_statusSym);
1659        slot_dup(ans, x, Matrix_pSym);
1660        slot_dup(ans, x, Matrix_xSym);
1661        INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1662        SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));
1663        REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;
1664        SET_SLOT(ans, Matrix_RXXSym, allocMatrix(REALSXP, 1, 1));
1665        REAL(GET_SLOT(ans, Matrix_RXXSym))[0] = NA_REAL;
1666        SET_SLOT(ans, Matrix_ZtXSym, allocMatrix(REALSXP, nz, 1));
1667        SET_SLOT(ans, Matrix_RZXSym, allocMatrix(REALSXP, nz, 1));
1668        LOGICAL(GET_SLOT(ans, Matrix_statusSym))[0] = 0;
1669        UNPROTECT(1);
1670        return ans;
1671      }      }
1672      UNPROTECT(2);
1673      return val;
1674    /**
1675     * Create an lme object from its components.  This is not done by
1676     * new("lme", ...) at the R level because of the possibility of
1677     * causing the copying of very large objects.
1678     *
1679     * @param call Pointer to the original call
1680     * @param facs pointer to the list of grouping factors
1681     * @param x pointer to the model matrices (may be of length zero)
1682     * @param model pointer to the model frame
1683     * @param REML pointer to a logical scalar indicating if REML is used
1684     * @param rep pointer to the converged ssclme object
1685     * @param fitted pointer to the fitted values
1686     * @param residuals pointer to the residuals
1687     *
1688     * @return an lme object
1689     */
1690    SEXP ssclme_to_lme(SEXP call, SEXP facs, SEXP x, SEXP model, SEXP REML,
1691                       SEXP rep, SEXP fitted, SEXP residuals)
1692    {
1693        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));
1694
1695        SET_SLOT(ans, install("call"), call);
1696        SET_SLOT(ans, install("facs"), facs);
1697        SET_SLOT(ans, Matrix_xSym, x);
1698        SET_SLOT(ans, install("model"), model);
1699        SET_SLOT(ans, install("REML"), REML);
1700        SET_SLOT(ans, install("rep"), rep);
1701        SET_SLOT(ans, install("fitted"), fitted);
1702        SET_SLOT(ans, install("residuals"), residuals);
1703        UNPROTECT(1);
1704        return ans;
1705  }  }
1706

Legend:
 Removed from v.100 changed lines Added in v.386