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 372, Sun Dec 5 13:53:29 2004 UTC
# 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.;  
   
     p = pp1 - 1;  
     b = RZX + p * n;  
     if (verbose) {  
         SEXP coef = PROTECT(ssclme_coef(x));  
         int lc = length(coef); double *cc = REAL(coef);  
1163    
         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        SEXP ans = PROTECT(allocMatrix(REALSXP, ncoef, ncoef)),
1513            base = PROTECT(ssclme_coef(x, Uncp)),
1514            current = PROTECT(duplicate(base)),
1515            gradient;
1516    
1517        for (j = 0; j < ncoef; j++) {
1518            double delta = (REAL(base)[j] ? 1.e-7 * REAL(base)[j] : 1.e-7);
1519            int i;
1520    
1521            for (i = 0; i < ncoef; i++) REAL(current)[i] = REAL(base)[i];
1522            REAL(current)[j] += delta/2.;
1523            ssclme_coefGets(x, current, Uncp);
1524            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1525            for (i = 0; i < ncoef; i++) REAL(ans)[j * ncoef + i] = REAL(gradient)[i];
1526            UNPROTECT(1);
1527            REAL(current)[j] -= delta;
1528            ssclme_coefGets(x, current, Uncp);
1529            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1530            for (i = 0; i < ncoef; i++)
1531                REAL(ans)[j * ncoef + i] =
1532                    (REAL(ans)[j * ncoef + i] - REAL(gradient)[i])/ delta;
1533            UNPROTECT(1);
1534            for (i = 0; i < j; i++) { /* symmetrize */
1535                REAL(ans)[j * ncoef + i] = REAL(ans)[i * ncoef + j] =
1536                    (REAL(ans)[j * ncoef + i] + REAL(ans)[i * ncoef + j])/2.;
1537            }
1538        }
1539        UNPROTECT(3);
1540        return ans;
1541    }
1542    
1543    /**
1544   * Calculate and return the fitted values.   * Calculate and return the fitted values.
1545   *   *
1546   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
# Line 1350  Line 1581 
1581              for (j = 0; j < nobs; ) {              for (j = 0; j < nobs; ) {
1582                  int nn = 1, lev = ff[j];                  int nn = 1, lev = ff[j];
1583                  /* check for adjacent rows with same factor level */                  /* check for adjacent rows with same factor level */
1584                  while (ff[j + nn] == lev) nn++;                  while ((j + nn) < nobs && ff[j + nn] == lev) nn++;
1585                  F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,                  F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
1586                                  &one, mm + j, &nobs,                                  &one, mm + j, &nobs,
1587                                  bb + (lev - 1) * nci, &nci,                                  bb + (lev - 1) * nci, &nci,
# Line 1471  Line 1702 
1702      UNPROTECT(1);      UNPROTECT(1);
1703      return ans;      return ans;
1704  }  }
1705    

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

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