SCM

SCM Repository

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

Diff of /pkg/src/ssclme.c

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

revision 176, Fri May 28 13:23:44 2004 UTC revision 565, Tue Feb 22 01:15:45 2005 UTC
# Line 97  Line 97 
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];
# Line 109  Line 110 
110                      }                      }
111                  }                  }
112              }              }
         }  
113          Free(map); Free(ncc);          Free(map); Free(ncc);
114      }      }
115  }  }
# Line 196  Line 196 
196      Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));      Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));
197      SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));      SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));
198      SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));      SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));
199      ldl_symbolic(nzcol, Ap, Ai, Lp, Parent,      R_ldl_symbolic(nzcol, Ap, Ai, Lp, Parent,
                  (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */  
                  (int *) R_alloc(nzcol, sizeof(int)), /* Flag */  
200                   (int *) NULL, (int *) NULL); /* P & Pinv */                   (int *) NULL, (int *) NULL); /* P & Pinv */
201      ssclme_calc_maxod(nzcol, Parent);      ssclme_calc_maxod(nzcol, Parent);
202      Lnz = Lp[nzcol];      Lnz = Lp[nzcol];
# Line 391  Line 389 
389                  nck = nc[k],                  nck = nc[k],
390                  scalar = ncj == 1 && nck == 1;                  scalar = ncj == 1 && nck == 1;
391              double              double
392                  *Zk = Z[k], *work;                  *Zk = Z[k], *work = (double *) NULL;
393    
394              if (!scalar) work = Calloc(ncj * nck, double);              if (!scalar) work = Calloc(ncj * nck, double);
395              for (i = 0; i < nobs; i++) {              for (i = 0; i < nobs; i++) {
396                  int ii, ind = -1, fpji = fpj[i] - 1,                  int ii, ind = -1, fpji = fpj[i] - 1,
# Line 477  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 567  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 625  Line 623 
623    
624   */   */
625  static  static
626  SEXP ldl_inverse(SEXP x)  void ldl_inverse(SEXP x)
627  {  {
628      SEXP      SEXP
629          Gpsl = GET_SLOT(x, Matrix_GpSym),          Gpsl = GET_SLOT(x, Matrix_GpSym),
# Line 705  Line 703 
703              Free(tmp); Free(ind);              Free(tmp); Free(ind);
704          }          }
705      }      }
     return R_NilValue;  
706  }  }
707    
708  /**  /**
# Line 747  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 831  Line 828 
828   *   *
829   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
830   *   *
831   * @return a vector containing the conditional modes of the random effects   * @return a list containing the conditional modes of the random effects
832   */   */
833  SEXP ssclme_ranef(SEXP x)  SEXP ssclme_ranef(SEXP x)
834  {  {
# Line 911  Line 908 
908   * order.   * order.
909   *   *
910   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
911     * @param Unconstr pointer to a logical object indicating if the
912     *        unconstrained parameterization should be used.
913   *   *
914   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of the values in the upper triangles of the
915   * Omega matrices   * Omega matrices
916   */   */
917  SEXP ssclme_coef(SEXP x)  SEXP ssclme_coef(SEXP x, SEXP Unconstr)
918  {  {
919      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
920      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
921          i, nf = length(Omega), vind;          i, nf = length(Omega), unc = asLogical(Unconstr), vind;
922      SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));      SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
923      double *vv = REAL(val);      double *vv = REAL(val);
924    
925      vind = 0;      vind = 0;                   /* index in vv */
926      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
927          int nci = nc[i];          int nci = nc[i], ncip1 = nci + 1;
928          if (nci == 1) {          if (nci == 1) {
929              vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];              vv[vind++] = (unc ?
930                              log(REAL(VECTOR_ELT(Omega, i))[0]) :
931                              REAL(VECTOR_ELT(Omega, i))[0]);
932          } else {          } else {
933              int j, k, odind = vind + nci, ncip1 = nci + 1;              if (unc) {          /* L log(D) L' factor of Omega[i,,] */
934                    int j, k, ncisq = nci * nci;
935                    double *tmp = Memcpy(Calloc(ncisq, double),
936                                         REAL(VECTOR_ELT(Omega, i)), ncisq);
937                    F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
938                    if (j)          /* should never happen */
939                        error("DPOTRF returned error code %d on Omega[[%d]]",
940                              j, i+1);
941                    for (j = 0; j < nci; j++) {
942                        double diagj = tmp[j * ncip1];
943                        vv[vind++] = 2. * log(diagj);
944                        for (k = j + 1; k < nci; k++) {
945                            tmp[j + k * nci] /= diagj;
946                        }
947                    }
948                    for (j = 0; j < nci; j++) {
949                        for (k = j + 1; k < nci; k++) {
950                            vv[vind++] = tmp[j + k * nci];
951                        }
952                    }
953                    Free(tmp);
954                } else {            /* upper triangle of Omega[i,,] */
955                    int j, k, odind = vind + nci;
956              double *omgi = REAL(VECTOR_ELT(Omega, i));              double *omgi = REAL(VECTOR_ELT(Omega, i));
957    
958              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
# Line 941  Line 964 
964              vind = odind;              vind = odind;
965          }          }
966      }      }
967        }
968      UNPROTECT(1);      UNPROTECT(1);
969      return val;      return val;
970  }  }
# Line 1055  Line 1079 
1079   *   *
1080   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1081   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
1082     * @param Unc pointer to a logical object indicating if the
1083     *        unconstrained parameterization should be used.
1084   *   *
1085   * @return R_NilValue   * @return R_NilValue
1086   */   */
1087  SEXP ssclme_coefGets(SEXP x, SEXP coef)  SEXP ssclme_coefGets(SEXP x, SEXP coef, SEXP Unc)
1088  {  {
1089      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1090      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1091            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1092          cind, i, nf = length(Omega),          cind, i, nf = length(Omega),
1093          *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));          unc = asLogical(Unc);
1094      double *cc = REAL(coef);      double *cc = REAL(coef);
1095    
1096      if (length(coef) != coef_length(nf, nc) || !isReal(coef))      if (length(coef) != coef_length(nf, nc) || !isReal(coef))
# Line 1073  Line 1100 
1100      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1101          int nci = nc[i];          int nci = nc[i];
1102          if (nci == 1) {          if (nci == 1) {
1103              REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];              REAL(VECTOR_ELT(Omega, i))[0] = (unc ?
1104                                                 exp(cc[cind++]) :
1105                                                 cc[cind++]);
1106          } else {          } else {
1107              int j, k, odind = cind + nci, ncip1 = nci + 1;              int odind = cind + nci, /* off-diagonal index */
1108              double *omgi = REAL(VECTOR_ELT(Omega, i));                  j, k,
1109                    ncip1 = nci + 1,
1110                    ncisq = nci * nci;
1111                double
1112                    *omgi = REAL(VECTOR_ELT(Omega, i));
1113                if (unc) {
1114                    double
1115                        *tmp = Calloc(ncisq, double),
1116                        diagj, one = 1., zero = 0.;
1117    
1118                    memset(omgi, 0, sizeof(double) * ncisq);
1119                    for (j = 0; j < nci; j++) {
1120                        tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1121                        for (k = j + 1; k < nci; k++) {
1122                            tmp[k*nci + j] = cc[odind++] * diagj;
1123                        }
1124                    }
1125                    F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1126                                    tmp, &nci, &zero, omgi, &nci);
1127                    Free(tmp);
1128                } else {
1129              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1130                  omgi[j * ncip1] = cc[cind++];                  omgi[j * ncip1] = cc[cind++];
1131                  for (k = j + 1; k < nci; k++) {                  for (k = j + 1; k < nci; k++) {
1132                      omgi[k*nci + j] = cc[odind++];                      omgi[k*nci + j] = cc[odind++];
1133                  }                  }
1134              }              }
1135                }
1136              cind = odind;              cind = odind;
1137          }          }
1138      }      }
# Line 1091  Line 1140 
1140      return x;      return x;
1141  }  }
1142    
1143    
1144  /**  /**
1145   * Perform a number of ECME steps for the REML or ML criterion.   * Returns the inverse of the updated Omega matrices for an ECME
1146     * iteration.  These matrices are also used in the gradient calculation.
1147   *   *
1148   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1149   * @param nsteps pointer to an integer scalar giving the number of ECME steps to perform   * @param REML indicator of REML being used
1150   * @param REMLp pointer to a logical scalar indicating if REML is to be used   * @param val pointer to a list of symmetric q_i by q_i matrices
  * @param verb pointer to a logical scalar indicating verbose mode  
  *  
  * @return NULL  
1151   */   */
1152  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)  static void
1153    common_ECME_gradient(SEXP x, int REML, SEXP val)
1154  {  {
1155      SEXP      SEXP
         Omega = GET_SLOT(x, Matrix_OmegaSym),  
1156          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ncsl = GET_SLOT(x, Matrix_ncSym),  
1157          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1158      int      int
1159          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1160          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1161          *nc = INTEGER(ncsl),          i, j, n = *INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1162          *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),          nf = length(val), nobs = nc[nf + 1], p = nc[nf] - 1;
         REML = asLogical(REMLp),  
         i, info, iter,  
         n = dims[0],  
         nEM = asInteger(nsteps),  
         nf = length(ncsl) - 2,  
         nobs = nc[nf + 1],  
         p,  
         pp1 = dims[1],  
         verbose = asLogical(verb);  
1163      double      double
1164          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
1165          *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),          *b = RZX + p * INTEGER(getAttrib(RZXsl, R_DimSymbol))[0],
1166          *b,          one = 1.0, zero = 0.0;
         alpha,  
         one = 1.,  
         zero = 0.;  
1167    
     p = pp1 - 1;  
     b = RZX + p * n;  
     if (verbose) {  
         SEXP coef = PROTECT(ssclme_coef(x));  
         int lc = length(coef); double *cc = REAL(coef);  
   
         ssclme_factor(x);  
         Rprintf("  EM iterations\n");  
         Rprintf("%3d %.3f", 0, dev[REML ? 1 : 0]);  
         for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);  
         Rprintf("\n");  
         UNPROTECT(1);  
     }  
     for (iter = 0; iter < nEM; iter++) {  
1168          ssclme_invert(x);          ssclme_invert(x);
1169          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
1170              int ki = Gp[i+1] - Gp[i],              int ki = Gp[i+1] - Gp[i],
1171                  nci = nc[i],                  nci = nc[i],
1172                  mi = ki/nci;                  mi = ki/nci;
1173              double              double
1174                  *vali = REAL(VECTOR_ELT(Omega, i));              *vali = REAL(VECTOR_ELT(val, i)),
1175                alpha = (double)(REML ? (nobs-p) : nobs)/((double) mi);
1176    
             alpha = ((double)(REML?(nobs-p):nobs))/((double)mi);  
1177              F77_CALL(dsyrk)("U", "N", &nci, &mi,              F77_CALL(dsyrk)("U", "N", &nci, &mi,
1178                              &alpha, b + Gp[i], &nci,                              &alpha, b + Gp[i], &nci,
1179                              &zero, vali, &nci);                              &zero, vali, &nci);
# Line 1161  Line 1182 
1182                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1183                              &one, vali, &nci);                              &one, vali, &nci);
1184              if (REML) {              if (REML) {
                 int j;  
1185                  for (j = 0; j < p; j++) {                  for (j = 0; j < p; j++) {
1186                      F77_CALL(dsyrk)("U", "N", &nci, &mi,                      F77_CALL(dsyrk)("U", "N", &nci, &mi,
1187                                  &alpha, RZX + Gp[i] + j*n, &nci,                                  &alpha, RZX + Gp[i] + j*n, &nci,
1188                                  &one, vali, &nci);                                  &one, vali, &nci);
1189                  }                  }
1190              }              }
             F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);  
             if (info)  
                 error("DPOTRF returned error code %d in Omega[[%d]] update",  
                       info, i + 1);  
             F77_CALL(dpotri)("U", &nci, vali, &nci, &info);  
             if (info)  
                 error("DPOTRI returned error code %d in Omega[[%d]] update",  
                       info, i + 1);  
1191          }          }
1192          status[0] = status[1] = 0;  }
1193          if (verbose) {  
1194              SEXP coef = PROTECT(ssclme_coef(x));  /**
1195              int lc = length(coef); double *cc = REAL(coef);   * Print the verbose output in the ECME iterations
1196     *
1197     * @param x pointer to an ssclme object
1198     * @param iter iteration number
1199     * @param REML indicator of whether REML is being used
1200     */
1201    static void EMsteps_verbose_print(SEXP x, int iter, int REML)
1202    {
1203        SEXP coef = PROTECT(ssclme_coef(x, ScalarLogical(0)));
1204        int i, lc = length(coef);
1205        double *cc = REAL(coef), *dev = REAL(GET_SLOT(x, Matrix_devianceSym));
1206    
1207    
1208              ssclme_factor(x);              ssclme_factor(x);
1209              Rprintf("%3d %.3f", iter + 1, dev[REML ? 1 : 0]);      if (iter == 0) Rprintf("  EM iterations\n");
1210        Rprintf("%3d %.3f", iter, dev[REML ? 1 : 0]);
1211              for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);              for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1212              Rprintf("\n");              Rprintf("\n");
1213              UNPROTECT(1);              UNPROTECT(1);
1214          }          }
1215    
1216    /**
1217     * Perform ECME steps for the REML or ML criterion.
1218     *
1219     * @param x pointer to an ssclme object
1220     * @param nsteps pointer to an integer scalar - the number of ECME steps to perform
1221     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1222     * @param verb pointer to a logical scalar indicating verbose mode
1223     *
1224     * @return NULL
1225     */
1226    SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1227    {
1228        SEXP
1229            Omega = GET_SLOT(x, Matrix_OmegaSym);
1230        int
1231            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1232            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1233            REML = asLogical(REMLp),
1234            i, info, iter,
1235            nEM = asInteger(nsteps),
1236            nf = length(Omega),
1237            verbose = asLogical(verb);
1238    
1239        if (verbose) EMsteps_verbose_print(x, 0, REML);
1240        for (iter = 0; iter < nEM; iter++) {
1241            common_ECME_gradient(x, REML, Omega);
1242            status[0] = status[1] = 0;
1243            for (i = 0; i < nf; i++) {
1244                double *vali = REAL(VECTOR_ELT(Omega, i));
1245    
1246                F77_CALL(dpotrf)("U", nc + i, vali, nc + i, &info);
1247                if (info)
1248                    error("DPOTRF returned error code %d in Omega[[%d]] update",
1249                          info, i + 1);
1250                F77_CALL(dpotri)("U", nc + i, vali, nc + i, &info);
1251                if (info)
1252                    error("DPOTRI returned error code %d in Omega[[%d]] update",
1253                          info, i + 1);
1254            }
1255            if (verbose) EMsteps_verbose_print(x, iter + 1, REML);
1256      }      }
1257      ssclme_factor(x);      ssclme_factor(x);
1258      return R_NilValue;      return R_NilValue;
1259  }  }
1260    
1261  /**  /**
1262     * Evaluate the gradient with respect to the indicators of the
1263     * positions in the Omega matrices.
1264     *
1265     * @param x Pointer to an ssclme object
1266     * @param REML indicator of whether REML is to be used
1267     * @param val Pointer to an object of the same structure as Omega
1268     */
1269    static void indicator_gradient(SEXP x, int REML, SEXP val)
1270    {
1271        SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1272        int
1273            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1274            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1275            i, nf = length(Omega);
1276    
1277        common_ECME_gradient(x, REML, val);
1278        for (i = 0; i < nf; i++) {
1279            int info, nci = nc[i], mi = (Gp[i+1] - Gp[i])/nci;
1280            double
1281                *work = Memcpy((double *) Calloc(nci * nci, double),
1282                               REAL(VECTOR_ELT(Omega, i)), nci * nci),
1283                alpha = (double) -mi, beta = (double) mi;
1284    
1285            F77_CALL(dpotrf)("U", &nci, work, &nci, &info);
1286            if (info)
1287                error("DPOTRF gave error code %d on Omega[[%d]]", info, i + 1);
1288            F77_CALL(dtrtri)("U", "N", &nci, work, &nci, &info);
1289            if (info)
1290                error("DPOTRI gave error code %d on Omega[[%d]]", info, i + 1);
1291            F77_CALL(dsyrk)("U", "N", &nci, &nci, &alpha, work, &nci,
1292                            &beta, REAL(VECTOR_ELT(val, i)), &nci);
1293            Free(work);
1294        }
1295    }
1296    
1297    /**
1298     * Overwrite a gradient with respect to positions in Omega[[i]] by the
1299     * gradient with respect to the unconstrained parameters.
1300     *
1301     * @param grad pointer to a gradient wrt positions.  Contents are overwritten.
1302     * @param Omega pointer to a list of symmetric matrices (upper storage).
1303     */
1304    static void unconstrained_gradient(SEXP grad, SEXP Omega)
1305    {
1306        int i, ii, j, nf = length(Omega);
1307        double one = 1.0, zero = 0.0;
1308    
1309        for (i = 0; i < nf; i++) {
1310            SEXP Omegai = VECTOR_ELT(Omega, i);
1311            int nci = *INTEGER(getAttrib(Omegai, R_DimSymbol)),
1312                ncisq = nci * nci, ncip1 = nci + 1;
1313            double *chol =
1314                Memcpy(Calloc(ncisq, double), REAL(Omegai), ncisq),
1315                *ai = REAL(VECTOR_ELT(grad, i)),
1316                *tmp = Calloc(ncisq, double);
1317    
1318            F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1319            if (j)
1320                error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1321            /* calculate and store grad[i,,] %*% t(R) */
1322            for (j = 0; j < nci; j++) {
1323                F77_CALL(dsymv)("U", &nci, &one, ai, &nci, chol + j, &nci,
1324                                &zero, tmp + j, &nci);
1325            }
1326            /* full symmetric product gives diagonals */
1327            F77_CALL(dtrmm)("R", "U", "T", "N", &nci, &nci, &one, chol, &nci,
1328                            Memcpy(ai, tmp, ncisq), &nci);
1329            /* overwrite upper triangle with gradients for positions in L' */
1330            for (ii = 1; ii < nci; ii++) {
1331                for (j = 0; j < ii; j++) {
1332                    ai[j + ii*nci] = 2. * chol[j*ncip1] * tmp[j + ii*nci];
1333                    ai[ii + j*nci] = 0.;
1334                }
1335            }
1336            Free(chol); Free(tmp);
1337        }
1338    }
1339    
1340    /**
1341     * Fills cvec with unlist(lapply(mList, function(el) el[upper.tri(el, strict = FALSE)]))
1342     *
1343     * @param mList pointer to a list of REAL matrices
1344     * @param nc number of columns in each matrix
1345     * @param cvec pointer to REAL vector to receive the result
1346     */
1347    static void upperTriList_to_vector(SEXP mList, int *nc, SEXP cvec)
1348    {
1349        int i, ii, j, nf = length(mList), pos = 0;
1350    
1351        pos = 0;                    /* position in vector */
1352        for (i = 0; i < nf; i++) {
1353            SEXP omgi = VECTOR_ELT(mList, i);
1354            int nci = nc[i];
1355            for (j = 0; j < nci; j++) {
1356                for (ii = 0; ii <= j; ii++) {
1357                    REAL(cvec)[pos++] = REAL(omgi)[ii + j * nci];
1358                }
1359            }
1360        }
1361    }
1362    
1363    /**
1364   * Return the gradient of the ML or REML deviance.   * Return the gradient of the ML or REML deviance.
1365   *   *
1366   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1367   * @param REMLp pointer to a logical scalar indicating if REML is to be used   * @param REMLp pointer to a logical scalar indicating if REML is to be used
1368   * @param Uncp pointer to a logical scalar indicating if the unconstrained parameterization is to be used   * @param Unc pointer to a logical scalar indicating if the
1369     *        unconstrained parameterization is to be used
1370     * @param OneVector pointer to a logical scalar indicating if result
1371     *        should be a single vector.
1372   *   *
1373   * @return pointer to a numeric vector of the gradient.   * @return pointer to the gradient as a list of matrices or as a vector.
1374   */   */
1375    SEXP ssclme_grad(SEXP x, SEXP REMLp, SEXP Unc, SEXP OneVector)
1376    {
1377        SEXP ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1378    
1379        indicator_gradient(x, asLogical(REMLp), ans);
1380        if (asLogical(Unc))
1381            unconstrained_gradient(ans, GET_SLOT(x, Matrix_OmegaSym));
1382        if (asLogical(OneVector)) {
1383            int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));
1384            SEXP val = PROTECT(allocVector(REALSXP, coef_length(length(ans), nc)));
1385    
1386            upperTriList_to_vector(ans, nc, val);
1387            UNPROTECT(2);
1388            return val;
1389        }
1390        UNPROTECT(1);
1391        return ans;
1392    }
1393    
1394  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1395  {  {
1396      SEXP      SEXP
# Line 1310  Line 1499 
1499  }  }
1500    
1501  /**  /**
1502     * Return the Hessian of the ML or REML deviance.  This is a
1503     * placeholder until I work out the evaluation of the analytic
1504     * Hessian, which probably will involve several helper functions.
1505     *
1506     * @param x pointer to an ssclme object
1507     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1508     * @param Uncp pointer to a logical scalar indicating if the
1509     * unconstrained parameterization is to be used
1510     *
1511     * @return pointer to an approximate Hessian matrix
1512     */
1513    SEXP ssclme_Hessian(SEXP x, SEXP REMLp, SEXP Uncp)
1514    {
1515        int j, ncoef = coef_length(length(GET_SLOT(x, Matrix_OmegaSym)),
1516                                   INTEGER(GET_SLOT(x, Matrix_ncSym)))
1517            /* , unc = asLogical(Uncp) */
1518            ;
1519        SEXP ans = PROTECT(allocMatrix(REALSXP, ncoef, ncoef)),
1520            base = PROTECT(ssclme_coef(x, Uncp)),
1521            current = PROTECT(duplicate(base)),
1522            gradient;
1523    
1524        for (j = 0; j < ncoef; j++) {
1525            double delta = (REAL(base)[j] ? 1.e-7 * REAL(base)[j] : 1.e-7);
1526            int i;
1527    
1528            for (i = 0; i < ncoef; i++) REAL(current)[i] = REAL(base)[i];
1529            REAL(current)[j] += delta/2.;
1530            ssclme_coefGets(x, current, Uncp);
1531            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1532            for (i = 0; i < ncoef; i++) REAL(ans)[j * ncoef + i] = REAL(gradient)[i];
1533            UNPROTECT(1);
1534            REAL(current)[j] -= delta;
1535            ssclme_coefGets(x, current, Uncp);
1536            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1537            for (i = 0; i < ncoef; i++)
1538                REAL(ans)[j * ncoef + i] =
1539                    (REAL(ans)[j * ncoef + i] - REAL(gradient)[i])/ delta;
1540            UNPROTECT(1);
1541            for (i = 0; i < j; i++) { /* symmetrize */
1542                REAL(ans)[j * ncoef + i] = REAL(ans)[i * ncoef + j] =
1543                    (REAL(ans)[j * ncoef + i] + REAL(ans)[i * ncoef + j])/2.;
1544            }
1545        }
1546        UNPROTECT(3);
1547        return ans;
1548    }
1549    
1550    /**
1551   * Calculate and return the fitted values.   * Calculate and return the fitted values.
1552   *   *
1553   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
# Line 1350  Line 1588 
1588              for (j = 0; j < nobs; ) {              for (j = 0; j < nobs; ) {
1589                  int nn = 1, lev = ff[j];                  int nn = 1, lev = ff[j];
1590                  /* check for adjacent rows with same factor level */                  /* check for adjacent rows with same factor level */
1591                  while (ff[j + nn] == lev) nn++;                  while ((j + nn) < nobs && ff[j + nn] == lev) nn++;
1592                  F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,                  F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
1593                                  &one, mm + j, &nobs,                                  &one, mm + j, &nobs,
1594                                  bb + (lev - 1) * nci, &nci,                                  bb + (lev - 1) * nci, &nci,
# Line 1452  Line 1690 
1690   * @param rep pointer to the converged ssclme object   * @param rep pointer to the converged ssclme object
1691   * @param fitted pointer to the fitted values   * @param fitted pointer to the fitted values
1692   * @param residuals pointer to the residuals   * @param residuals pointer to the residuals
1693     * @param terms pointer to a terms object (redundant if model is non-empty)
1694     * @param assign pointer to an integer assign vector
1695   *   *
1696   * @return an lme object   * @return an lme object
1697   */   */
1698  SEXP ssclme_to_lme(SEXP call, SEXP facs, SEXP x, SEXP model, SEXP REML,  SEXP ssclme_to_lme(SEXP call, SEXP facs, SEXP x, SEXP model, SEXP REML,
1699                     SEXP rep, SEXP fitted, SEXP residuals)                     SEXP rep, SEXP fitted, SEXP residuals, SEXP terms,
1700                       SEXP assign)
1701  {  {
1702      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));
1703    
# Line 1468  Line 1709 
1709      SET_SLOT(ans, install("rep"), rep);      SET_SLOT(ans, install("rep"), rep);
1710      SET_SLOT(ans, install("fitted"), fitted);      SET_SLOT(ans, install("fitted"), fitted);
1711      SET_SLOT(ans, install("residuals"), residuals);      SET_SLOT(ans, install("residuals"), residuals);
1712        SET_SLOT(ans, install("terms"), terms);
1713        SET_SLOT(ans, install("assign"), assign);
1714      UNPROTECT(1);      UNPROTECT(1);
1715      return ans;      return ans;
1716  }  }
1717    

Legend:
Removed from v.176  
changed lines
  Added in v.565

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