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 397, Thu Dec 16 19:15:32 2004 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 915  Line 912 
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 nci = nc[i];          int nci = nc[i], ncip1 = nci + 1;
926          if (nci == 1) {          if (nci == 1) {
927              vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];              vv[vind++] = (unc ?
928                              log(REAL(VECTOR_ELT(Omega, i))[0]) :
929                              REAL(VECTOR_ELT(Omega, i))[0]);
930          } else {          } else {
931              int j, k, odind = vind + nci, ncip1 = nci + 1;              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++) {
# Line 941  Line 962 
962              vind = odind;              vind = odind;
963          }          }
964      }      }
965        }
966      UNPROTECT(1);      UNPROTECT(1);
967      return val;      return val;
968  }  }
# Line 1058  Line 1080 
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 1073  Line 1096 
1096      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1097          int nci = nc[i];          int nci = nc[i];
1098          if (nci == 1) {          if (nci == 1) {
1099              REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];              REAL(VECTOR_ELT(Omega, i))[0] = (unc ?
1100                                                 exp(cc[cind++]) :
1101                                                 cc[cind++]);
1102          } else {          } else {
1103              int j, k, odind = cind + nci, ncip1 = nci + 1;              int odind = cind + nci, /* off-diagonal index */
1104              double *omgi = REAL(VECTOR_ELT(Omega, i));                  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++) {
1116                        tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1117                        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++) {              for (j = 0; j < nci; j++) {
1126                  omgi[j * ncip1] = cc[cind++];                  omgi[j * ncip1] = cc[cind++];
1127                  for (k = j + 1; k < nci; k++) {                  for (k = j + 1; k < nci; k++) {
1128                      omgi[k*nci + j] = cc[odind++];                      omgi[k*nci + j] = cc[odind++];
1129                  }                  }
1130              }              }
1131                }
1132              cind = odind;              cind = odind;
1133          }          }
1134      }      }
# Line 1091  Line 1136 
1136      return x;      return x;
1137  }  }
1138    
1139    
1140  /**  /**
1141   * Perform a number of ECME steps for the REML or ML criterion.   * 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   * @param x pointer to an ssclme object
1145   * @param nsteps pointer to an integer scalar giving the number of ECME steps to perform   * @param REML indicator of REML being used
1146   * @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  
1147   */   */
1148  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)  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) {  
         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++) {  
1164          ssclme_invert(x);          ssclme_invert(x);
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 1161  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 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);  
1187          }          }
1188          status[0] = status[1] = 0;  }
1189          if (verbose) {  
1190              SEXP coef = PROTECT(ssclme_coef(x));  /**
1191              int lc = length(coef); double *cc = REAL(coef);   * 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);              ssclme_factor(x);
1205              Rprintf("%3d %.3f", iter + 1, dev[REML ? 1 : 0]);      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]);              for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1208              Rprintf("\n");              Rprintf("\n");
1209              UNPROTECT(1);              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++) {
1237            common_ECME_gradient(x, REML, Omega);
1238            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  /**  /**
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    
1273        common_ECME_gradient(x, REML, val);
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     */
1300    static void unconstrained_gradient(SEXP grad, SEXP Omega)
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),
1311                *ai = REAL(VECTOR_ELT(grad, i)),
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.   * Return the gradient of the ML or REML deviance.
1361   *   *
1362   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1363   * @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
1364   * @param Uncp pointer to a logical scalar indicating if the unconstrained parameterization is to be used   * @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 a numeric vector of the gradient.   * @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    
1373        indicator_gradient(x, asLogical(REMLp), ans);
1374        if (asLogical(Unc))
1375            unconstrained_gradient(ans, GET_SLOT(x, Matrix_OmegaSym));
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)  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1389  {  {
1390      SEXP      SEXP
# Line 1310  Line 1493 
1493  }  }
1494    
1495  /**  /**
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)),
1516            gradient;
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);
1525            PROTECT(gradient = ssclme_gradient(x, REMLp, 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);
1530            PROTECT(gradient = ssclme_gradient(x, REMLp, 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.   * Calculate and return the fitted values.
1546   *   *
1547   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
# Line 1350  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 1452  Line 1684 
1684   * @param rep pointer to the converged ssclme object   * @param rep pointer to the converged ssclme object
1685   * @param fitted pointer to the fitted values   * @param fitted pointer to the fitted values
1686   * @param residuals pointer to the residuals   * @param residuals pointer to the residuals
1687     * @param terms pointer to a terms object (redundant if model is non-empty)
1688     * @param assign pointer to an integer assign vector
1689     * @param contrasts pointer to a list of contrast function names
1690   *   *
1691   * @return an lme object   * @return an lme object
1692   */   */
1693  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,
1694                     SEXP rep, SEXP fitted, SEXP residuals)                     SEXP rep, SEXP fitted, SEXP residuals, SEXP terms,
1695                       SEXP assign, SEXP contrasts)
1696  {  {
1697      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));
1698    
# Line 1468  Line 1704 
1704      SET_SLOT(ans, install("rep"), rep);      SET_SLOT(ans, install("rep"), rep);
1705      SET_SLOT(ans, install("fitted"), fitted);      SET_SLOT(ans, install("fitted"), fitted);
1706      SET_SLOT(ans, install("residuals"), residuals);      SET_SLOT(ans, install("residuals"), residuals);
1707        SET_SLOT(ans, install("terms"), terms);
1708        SET_SLOT(ans, install("assign"), assign);
1709        SET_SLOT(ans, install("contrasts"), contrasts);
1710      UNPROTECT(1);      UNPROTECT(1);
1711      return ans;      return ans;
1712  }  }
1713    

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

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