SCM

SCM Repository

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

Diff of /pkg/src/lmer.c

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

revision 752, Sun May 29 02:14:58 2005 UTC revision 753, Sun May 29 02:15:35 2005 UTC
# Line 1  Line 1 
1  #include "lmer.h"  #include "lmer.h"
 /* TODO  
  * - The egsingle example with ~year|childid+schoolid shows an unusual  
  *   drop in the deviance when switching from ECME to optim.  Is it real?  
  *   (Apparently so.)  
  * - Remove the fill_nnz function and the PAR and BLK macros.  Do the  
  *   allocation of temporary storage once only.  
  * - Perhaps change lmer_crosstab to omit to counts and work with  
  *   logical sparse matrix forms only.  (We only use the form, not the  
  *   counts)  
  */  
2    
3  /**  /**
4   * Check validity of an lmer object.   * Check validity of an lmer object.
# Line 199  Line 189 
189  void  void
190  lmer_populate(SEXP val)  lmer_populate(SEXP val)
191  {  {
192      SEXP D, L, Parent, ZZpO, flist = GET_SLOT(val, Matrix_flistSym),      SEXP D, L, Parent, ZZpO,
193            flist = GET_SLOT(val, Matrix_flistSym),
194          perm, Omega, ZtZ = GET_SLOT(val, Matrix_ZtZSym);          perm, Omega, ZtZ = GET_SLOT(val, Matrix_ZtZSym);
195      SEXP fnms = getAttrib(flist, R_NamesSymbol);      SEXP fnms = getAttrib(flist, R_NamesSymbol);
196      int j, k, nf = length(flist);      int j, k, nf = length(flist);
# Line 233  Line 224 
224          nlev[j] = length(getAttrib(VECTOR_ELT(flist, j), R_LevelsSymbol));          nlev[j] = length(getAttrib(VECTOR_ELT(flist, j), R_LevelsSymbol));
225          Gp[j + 1] = Gp[j] + nc[j] * nlev[j];          Gp[j + 1] = Gp[j] + nc[j] * nlev[j];
226          SET_VECTOR_ELT(D, j, alloc3Darray(REALSXP, nc[j], nc[j], nlev[j]));          SET_VECTOR_ELT(D, j, alloc3Darray(REALSXP, nc[j], nc[j], nlev[j]));
227            AZERO(REAL(VECTOR_ELT(D, j)), nc[j] * nc[j] * nlev[j]);
228          SET_VECTOR_ELT(Omega, j, allocMatrix(REALSXP, nc[j], nc[j]));          SET_VECTOR_ELT(Omega, j, allocMatrix(REALSXP, nc[j], nc[j]));
229          SET_VECTOR_ELT(ZZpO, j, duplicate(VECTOR_ELT(ZtZ, Lind(j, j))));          SET_VECTOR_ELT(ZZpO, j, duplicate(VECTOR_ELT(ZtZ, Lind(j, j))));
230          for (k = j; k < nf; k++)          for (k = j; k < nf; k++)
# Line 977  Line 969 
969    
970          alloc_tmp_ind(nf, nc, nlevs, ParP, tmp, ind);          alloc_tmp_ind(nf, nc, nlevs, ParP, tmp, ind);
971          /* Create bVar arrays as crossprod of column blocks of D^{-T/2}%*%L^{-1} */          /* Create bVar arrays as crossprod of column blocks of D^{-T/2}%*%L^{-1} */
972          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) { /* ith column of outer blocks */
973              int j, k, kj, nci = nc[i];              int j, k, kj, nci = nc[i];
974              int ncisqr = nci * nci;              int ncisqr = nci * nci;
975              double *Di = REAL(VECTOR_ELT(DP, i)),              double *Di = REAL(VECTOR_ELT(DP, i)),
976                  *bVi = REAL(VECTOR_ELT(bVarP, i));                  *bVi = REAL(VECTOR_ELT(bVarP, i));
977    
978              AZERO(bVi, nlevs[i] * ncisqr); /* zero the accumulator */              AZERO(bVi, ncisqr * nlevs[i]);
979              for (j = 0; j < nlevs[i]; j++) {              for (j = 0; j < nlevs[i]; j++) {
980                  double *bVij = bVi + j * ncisqr, *Dij = Di + j * ncisqr;                  double *bVij = bVi + j * ncisqr, *Dij = Di + j * ncisqr;
981    
# Line 991  Line 983 
983                                  &nci, &zero, bVij, &nci);                                  &nci, &zero, bVij, &nci);
984                  /* count non-zero blocks; allocate and zero storage */                  /* count non-zero blocks; allocate and zero storage */
985                  fill_ind(i, j, nf, ParP, nnz, ind);                  fill_ind(i, j, nf, ParP, nnz, ind);
986                  /* kth row of outer blocks */  
987                  for (k = i; k < nf; k++) {                  for (k = i; k < nf; k++) { /* kth row of outer blocks */
988                      SEXP Lki = VECTOR_ELT(LP, Lind(k, i));                      SEXP Lki = VECTOR_ELT(LP, Lind(k, i));
989                      int *Lkii = INTEGER(GET_SLOT(Lki, Matrix_iSym)),                      int *Lkii = INTEGER(GET_SLOT(Lki, Matrix_iSym)),
990                          *Lkip = INTEGER(GET_SLOT(Lki, Matrix_pSym));                          *Lkip = INTEGER(GET_SLOT(Lki, Matrix_pSym));
# Line 1015  Line 1007 
1007                          Lkip = INTEGER(GET_SLOT(Lki, Matrix_pSym));                          Lkip = INTEGER(GET_SLOT(Lki, Matrix_pSym));
1008                          Lkix = REAL(GET_SLOT(Lki, Matrix_xSym));                          Lkix = REAL(GET_SLOT(Lki, Matrix_xSym));
1009                          for (kj = 0; kj < nnz[kk]; kj++) {                          for (kj = 0; kj < nnz[kk]; kj++) {
1010                              int col = ind[kk][kj], k1;                              int col = ind[kk][kj], k1, szkk = nc[i] * nc[kk];
1011    
1012                              for (k1 = Lkip[col]; k1 < Lkip[col + 1]; k1++) {                              for (k1 = Lkip[col]; k1 < Lkip[col + 1]; k1++) {
1013                                  if ((kk == k) && col >= Lkii[k1]) break;                                  if ((kk == k) && col >= Lkii[k1]) break;
1014                                  F77_CALL(dgemm)("N", "N", nc+k, &nci, nc+kk,                                  F77_CALL(dgemm)("N", "N", &nc[k], &nci, &nc[kk],
1015                                                  &minus1, Lkix + k1 * szk,                                                  &minus1, Lkix + k1 * szk,
1016                                                  nc + k, tmp[kk] + kj * szk,                                                  &nc[k], tmp[kk] + kj * szkk,
1017                                                  nc + k, &one,                                                  &nc[kk], &one,
1018                                                  tmp[k]+fsrch(Lkii[k1],ind[k],nnz[k])*sz,                                                  tmp[k] +
1019                                                  nc + k);                                                  fsrch(Lkii[k1],ind[k],nnz[k])*sz,
1020                                                    &nc[k]);
1021                              }                              }
1022                          }                          }
1023                      }                      }
# Line 1094  Line 1087 
1087  }  }
1088    
1089  /**  /**
1090   * Extract the upper triangles of the Omega matrices.  These aren't   * Extract parameters from the Omega matrices.  These aren't
1091   * "coefficients" but the extractor is called coef for historical   * "coefficients" but the extractor is called coef for historical
1092   * reasons.  Within each group these values are in the order of the   * reasons.  Within each group these values are in the order of the
1093   * diagonal entries first then the strict upper triangle in row   * diagonal entries first then the strict upper triangle in row
1094   * order.   * order.
1095   *   *
1096     * The parameters can be returned in three forms:
1097     *   0 - nonlinearly constrained - elements of the relative precision matrix
1098     *   1 - unconstrained - from the LDL' decomposition - logarithms of
1099     *       the diagonal elements of D
1100     *   2 - box constrained - also from the LDL' decomposition - inverses
1101     *       of the diagonal elements of D
1102     *
1103   * @param x pointer to an lme object   * @param x pointer to an lme object
1104   * @param Unc pointer to a logical scalar indicating if the parameters   * @param pType pointer to an integer scalar indicating the form of the
1105   * are in the unconstrained form.   *        parameters to be returned.
1106   *   *
1107   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of the values in the upper triangles of the
1108   * Omega matrices   * Omega matrices
1109   */   */
1110  SEXP lmer_coef(SEXP x, SEXP Unc)  SEXP lmer_coef(SEXP x, SEXP pType)
1111  {  {
1112      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1113      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1114          i, nf = length(Omega), unc = asLogical(Unc), vind;          i, nf = length(Omega), ptyp = asInteger(pType), vind;
1115      SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));      SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1116      double *vv = REAL(val);      double *vv = REAL(val);
1117    
# Line 1119  Line 1119 
1119      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1120          int nci = nc[i], ncip1 = nci + 1;          int nci = nc[i], ncip1 = nci + 1;
1121          if (nci == 1) {          if (nci == 1) {
1122              vv[vind++] = (unc ?              double dd = REAL(VECTOR_ELT(Omega, i))[0];
1123                            log(REAL(VECTOR_ELT(Omega, i))[0]) :              vv[vind++] = ptyp ? ((ptyp == 1) ? log(dd) : 1./dd) : dd;
                           REAL(VECTOR_ELT(Omega, i))[0]);  
1124          } else {          } else {
1125              if (unc) {          /* L log(D) L' factor of Omega[,,i] */              if (ptyp) { /* L log(D) L' factor of Omega[,,i] */
1126                  int j, k, ncisq = nci * nci;                  int j, k, ncisq = nci * nci;
1127                  double *tmp = Memcpy(Calloc(ncisq, double),                  double *tmp = Memcpy(Calloc(ncisq, double),
1128                                       REAL(VECTOR_ELT(Omega, i)), ncisq);                                       REAL(VECTOR_ELT(Omega, i)), ncisq);
# Line 1133  Line 1132 
1132                            j, i+1);                            j, i+1);
1133                  for (j = 0; j < nci; j++) {                  for (j = 0; j < nci; j++) {
1134                      double diagj = tmp[j * ncip1];                      double diagj = tmp[j * ncip1];
1135                      vv[vind++] = 2. * log(diagj);                      vv[vind++] = (ptyp == 1) ? (2. * log(diagj)) :
1136                            1./(diagj * diagj);
1137                      for (k = j + 1; k < nci; k++) {                      for (k = j + 1; k < nci; k++) {
1138                          tmp[j + k * nci] /= diagj;                          tmp[j + k * nci] /= diagj;
1139                      }                      }
# Line 1164  Line 1164 
1164    
1165    
1166  /**  /**
1167   * Assign the upper triangles of the Omega matrices.   * Assign the upper triangles of the Omega matrices according to a
1168   * (Called coef for historical reasons.)   * vector of parameters.
1169   *   *
1170   * @param x pointer to an lme object   * @param x pointer to an lme object
1171   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
1172   * @param Unc pointer to a logical scalar indicating if the parameters   * @param pType pointer to an integer scalar
  * are in the unconstrained form.  
1173   *   *
1174   * @return R_NilValue   * @return R_NilValue
1175   */   */
1176  SEXP lmer_coefGets(SEXP x, SEXP coef, SEXP Unc)  SEXP lmer_coefGets(SEXP x, SEXP coef, SEXP pType)
1177  {  {
1178      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1179      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1180          *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),          *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1181          cind, i, nf = length(Omega),          cind, i, nf = length(Omega),
1182          unc = asLogical(Unc);          ptyp = asInteger(pType);
1183      double *cc = REAL(coef);      double *cc = REAL(coef);
1184    
1185      if (length(coef) != coef_length(nf, nc) || !isReal(coef))      if (length(coef) != coef_length(nf, nc) || !isReal(coef))
# Line 1190  Line 1189 
1189      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1190          int nci = nc[i];          int nci = nc[i];
1191          if (nci == 1) {          if (nci == 1) {
1192              REAL(VECTOR_ELT(Omega, i))[0] = (unc ?              double dd = cc[cind++];
1193                                               exp(cc[cind++]) :              REAL(VECTOR_ELT(Omega, i))[0] =
1194                                               cc[cind++]);                  ptyp ? ((ptyp == 1) ? exp(dd) : 1./dd) : dd;
1195          } else {          } else {
1196              int odind = cind + nci, /* off-diagonal index */              int odind = cind + nci, /* off-diagonal index */
1197                  j, k,                  j, k,
# Line 1200  Line 1199 
1199                  ncisq = nci * nci;                  ncisq = nci * nci;
1200              double              double
1201                  *omgi = REAL(VECTOR_ELT(Omega, i));                  *omgi = REAL(VECTOR_ELT(Omega, i));
1202              if (unc) {              if (ptyp) {
1203                  double                  double *tmp = Calloc(ncisq, double),
                     *tmp = Calloc(ncisq, double),  
1204                      diagj, one = 1., zero = 0.;                      diagj, one = 1., zero = 0.;
1205    
1206                  AZERO(omgi, ncisq);                  AZERO(omgi, ncisq);
1207                  for (j = 0; j < nci; j++) {                  for (j = 0; j < nci; j++) {
1208                      tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);                      double dd = cc[cind++];
1209                        tmp[j * ncip1] = diagj =
1210                            (ptyp == 1) ? exp(dd/2.) : sqrt(1./dd);
1211                      for (k = j + 1; k < nci; k++) {                      for (k = j + 1; k < nci; k++) {
1212                          tmp[k*nci + j] = cc[odind++] * diagj;                          tmp[k*nci + j] = cc[odind++] * diagj;
1213                      }                      }
# Line 1560  Line 1560 
1560      return val;      return val;
1561  }  }
1562    
1563  SEXP lmer_gradient(SEXP x, SEXP REMLp, SEXP Uncp)  SEXP lmer_gradient(SEXP x, SEXP REMLp, SEXP pType)
1564  {  {
1565      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1566      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1567          dind, i, ifour = 4, info, ione = 1, nf = length(Omega),          dind, i, ifour = 4, info, ione = 1, nf = length(Omega),
1568          odind, unc = asLogical(Uncp);          odind, ptyp = asInteger(pType);
1569      SEXP      SEXP
1570          firstDer = lmer_firstDer(x, PROTECT(EM_grad_array(nf, nc))),          firstDer = lmer_firstDer(x, PROTECT(EM_grad_array(nf, nc))),
1571          val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));          val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
# Line 1585  Line 1585 
1585                          REAL(VECTOR_ELT(firstDer, i)), &ncisqr,                          REAL(VECTOR_ELT(firstDer, i)), &ncisqr,
1586                          cc, &ifour, &zero, tmp, &ncisqr);                          cc, &ifour, &zero, tmp, &ncisqr);
1587          if (nci == 1) {          if (nci == 1) {
1588              REAL(val)[dind++] = (unc ? Omgi[0] : 1.) * tmp[0];              REAL(val)[dind++] =
1589                    (ptyp?((ptyp == 1)?Omgi[0]: -Omgi[0] * Omgi[0]) : 1) * tmp[0];
1590          } else {          } else {
1591              int ii, j, ncip1 = nci + 1;              int ii, j, ncip1 = nci + 1;
1592    
1593              odind = dind + nci; /* index into val for off-diagonals */              odind = dind + nci; /* index into val for off-diagonals */
1594              if (unc) {              if (ptyp) {
1595                  double *chol = Memcpy(Calloc(ncisqr, double),                  double *chol = Memcpy(Calloc(ncisqr, double),
1596                                        REAL(VECTOR_ELT(Omega, i)), ncisqr),                                        REAL(VECTOR_ELT(Omega, i)), ncisqr),
1597                      *tmp2 = Calloc(ncisqr, double);                      *tmp2 = Calloc(ncisqr, double);
# Line 1615  Line 1616 
1616                          tmp[ii + j*nci] = 0.;                          tmp[ii + j*nci] = 0.;
1617                      }                      }
1618                  }                  }
1619                  Free(chol); Free(tmp2);                  if (ptyp > 1)
1620                        for (ii = 0; ii < nci; ii++) {
1621                            int ind = ii * ncip1;
1622                            double sqrtd = chol[ind];
1623                            tmp[ind] *= -(sqrtd*sqrtd);
1624                        }
1625              }              }
1626              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1627                  REAL(val)[dind + j] = tmp[j * ncip1];                  REAL(val)[dind + j] = tmp[j * ncip1];

Legend:
Removed from v.752  
changed lines
  Added in v.753

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