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 247, Mon Jul 12 17:41:57 2004 UTC
# Line 391  Line 391 
391                  nck = nc[k],                  nck = nc[k],
392                  scalar = ncj == 1 && nck == 1;                  scalar = ncj == 1 && nck == 1;
393              double              double
394                  *Zk = Z[k], *work;                  *Zk = Z[k], *work = (double *) NULL;
395    
396              if (!scalar) work = Calloc(ncj * nck, double);              if (!scalar) work = Calloc(ncj * nck, double);
397              for (i = 0; i < nobs; i++) {              for (i = 0; i < nobs; i++) {
398                  int ii, ind = -1, fpji = fpj[i] - 1,                  int ii, ind = -1, fpji = fpj[i] - 1,
# Line 625  Line 626 
626    
627   */   */
628  static  static
629  SEXP ldl_inverse(SEXP x)  void ldl_inverse(SEXP x)
630  {  {
631      SEXP      SEXP
632          Gpsl = GET_SLOT(x, Matrix_GpSym),          Gpsl = GET_SLOT(x, Matrix_GpSym),
# Line 705  Line 706 
706              Free(tmp); Free(ind);              Free(tmp); Free(ind);
707          }          }
708      }      }
     return R_NilValue;  
709  }  }
710    
711  /**  /**
# Line 915  Line 915 
915   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of the values in the upper triangles of the
916   * Omega matrices   * Omega matrices
917   */   */
918  SEXP ssclme_coef(SEXP x)  SEXP ssclme_coef(SEXP x, SEXP Unconstr)
919  {  {
920      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
921      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
922          i, nf = length(Omega), vind;          i, nf = length(Omega), unc = asLogical(Unconstr), vind;
923      SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));      SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
924      double *vv = REAL(val);      double *vv = REAL(val);
925    
926      vind = 0;      vind = 0;                   /* index in vv */
927      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
928          int nci = nc[i];          int nci = nc[i], ncip1 = nci + 1;
929          if (nci == 1) {          if (nci == 1) {
930              vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];              vv[vind++] = (unc ?
931                              log(REAL(VECTOR_ELT(Omega, i))[0]) :
932                              REAL(VECTOR_ELT(Omega, i))[0]);
933          } else {          } else {
934              int j, k, odind = vind + nci, ncip1 = nci + 1;              if (unc) {          /* L log(D) L' factor of Omega[i,,] */
935                    int j, k, ncisq = nci * nci;
936                    double *tmp = Memcpy(Calloc(ncisq, double),
937                                         REAL(VECTOR_ELT(Omega, i)), ncisq);
938                    F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
939                    if (j)          /* should never happen */
940                        error("DPOTRF returned error code %d on Omega[[%d]]",
941                              j, i+1);
942                    for (j = 0; j < nci; j++) {
943                        double diagj = tmp[j * ncip1];
944                        vv[vind++] = 2. * log(diagj);
945                        for (k = j + 1; k < nci; k++) {
946                            tmp[j + k * nci] /= diagj;
947                        }
948                    }
949                    for (j = 0; j < nci; j++) {
950                        for (k = j + 1; k < nci; k++) {
951                            vv[vind++] = tmp[j + k * nci];
952                        }
953                    }
954                    Free(tmp);
955                } else {            /* upper triangle of Omega[i,,] */
956                    int j, k, odind = vind + nci;
957              double *omgi = REAL(VECTOR_ELT(Omega, i));              double *omgi = REAL(VECTOR_ELT(Omega, i));
958    
959              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
# Line 941  Line 965 
965              vind = odind;              vind = odind;
966          }          }
967      }      }
968        }
969      UNPROTECT(1);      UNPROTECT(1);
970      return val;      return val;
971  }  }
# Line 1058  Line 1083 
1083   *   *
1084   * @return R_NilValue   * @return R_NilValue
1085   */   */
1086  SEXP ssclme_coefGets(SEXP x, SEXP coef)  SEXP ssclme_coefGets(SEXP x, SEXP coef, SEXP Unc)
1087  {  {
1088      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1089      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1090            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1091          cind, i, nf = length(Omega),          cind, i, nf = length(Omega),
1092          *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));          unc = asLogical(Unc);
1093      double *cc = REAL(coef);      double *cc = REAL(coef);
1094    
1095      if (length(coef) != coef_length(nf, nc) || !isReal(coef))      if (length(coef) != coef_length(nf, nc) || !isReal(coef))
# Line 1073  Line 1099 
1099      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1100          int nci = nc[i];          int nci = nc[i];
1101          if (nci == 1) {          if (nci == 1) {
1102              REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];              REAL(VECTOR_ELT(Omega, i))[0] = (unc ?
1103                                                 exp(cc[cind++]) :
1104                                                 cc[cind++]);
1105          } else {          } else {
1106              int j, k, odind = cind + nci, ncip1 = nci + 1;              int odind = cind + nci, /* off-diagonal index */
1107              double *omgi = REAL(VECTOR_ELT(Omega, i));                  j, k,
1108                    ncip1 = nci + 1,
1109                    ncisq = nci * nci;
1110                double
1111                    *omgi = REAL(VECTOR_ELT(Omega, i));
1112                if (unc) {
1113                    double
1114                        *tmp = Calloc(ncisq, double),
1115                        diagj, one = 1., zero = 0.;
1116    
1117                    memset(omgi, 0, sizeof(double) * ncisq);
1118                    for (j = 0; j < nci; j++) {
1119                        tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1120                        for (k = j + 1; k < nci; k++) {
1121                            tmp[k*nci + j] = cc[odind++] * diagj;
1122                        }
1123                    }
1124                    F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1125                                    tmp, &nci, &zero, omgi, &nci);
1126                    Free(tmp);
1127                } else {
1128              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1129                  omgi[j * ncip1] = cc[cind++];                  omgi[j * ncip1] = cc[cind++];
1130                  for (k = j + 1; k < nci; k++) {                  for (k = j + 1; k < nci; k++) {
1131                      omgi[k*nci + j] = cc[odind++];                      omgi[k*nci + j] = cc[odind++];
1132                  }                  }
1133              }              }
1134                }
1135              cind = odind;              cind = odind;
1136          }          }
1137      }      }
# Line 1091  Line 1139 
1139      return x;      return x;
1140  }  }
1141    
1142    
1143  /**  /**
1144   * Perform a number of ECME steps for the REML or ML criterion.   * Returns the inverse of the updated Omega matrices for an ECME
1145     * iteration.  These matrices are also used in the gradient calculation.
1146   *   *
1147   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1148   * @param nsteps pointer to an integer scalar giving the number of ECME steps to perform   * @param REML indicator of REML being used
1149   * @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  
1150   */   */
1151  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)  static void
1152    common_ECME_gradient(SEXP x, int REML, SEXP val)
1153  {  {
1154      SEXP      SEXP
         Omega = GET_SLOT(x, Matrix_OmegaSym),  
1155          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ncsl = GET_SLOT(x, Matrix_ncSym),  
1156          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1157      int      int
1158          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1159          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1160          *nc = INTEGER(ncsl),          i, j, n = *INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1161          *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);  
1162      double      double
1163          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
1164          *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),          *b = RZX + p * INTEGER(getAttrib(RZXsl, R_DimSymbol))[0],
1165          *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);  
1166    
         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++) {  
1167          ssclme_invert(x);          ssclme_invert(x);
1168          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
1169              int ki = Gp[i+1] - Gp[i],              int ki = Gp[i+1] - Gp[i],
1170                  nci = nc[i],                  nci = nc[i],
1171                  mi = ki/nci;                  mi = ki/nci;
1172              double              double
1173                  *vali = REAL(VECTOR_ELT(Omega, i));              *vali = REAL(VECTOR_ELT(val, i)),
1174                alpha = (double)(REML ? (nobs-p) : nobs)/((double) mi);
1175    
             alpha = ((double)(REML?(nobs-p):nobs))/((double)mi);  
1176              F77_CALL(dsyrk)("U", "N", &nci, &mi,              F77_CALL(dsyrk)("U", "N", &nci, &mi,
1177                              &alpha, b + Gp[i], &nci,                              &alpha, b + Gp[i], &nci,
1178                              &zero, vali, &nci);                              &zero, vali, &nci);
# Line 1161  Line 1181 
1181                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1182                              &one, vali, &nci);                              &one, vali, &nci);
1183              if (REML) {              if (REML) {
                 int j;  
1184                  for (j = 0; j < p; j++) {                  for (j = 0; j < p; j++) {
1185                      F77_CALL(dsyrk)("U", "N", &nci, &mi,                      F77_CALL(dsyrk)("U", "N", &nci, &mi,
1186                                  &alpha, RZX + Gp[i] + j*n, &nci,                                  &alpha, RZX + Gp[i] + j*n, &nci,
1187                                  &one, vali, &nci);                                  &one, vali, &nci);
1188                  }                  }
1189              }              }
             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);  
1190          }          }
1191          status[0] = status[1] = 0;  }
1192          if (verbose) {  
1193              SEXP coef = PROTECT(ssclme_coef(x));  /**
1194     * Print the verbose output in the ECME iterations
1195     *
1196     * @param x pointer to an ssclme object
1197     * @param iter iteration number
1198     * @param REML indicator of whether REML is being used
1199     */
1200    static void EMsteps_verbose_print(SEXP x, int iter, int REML)
1201    {
1202        SEXP coef = PROTECT(ssclme_coef(x, ScalarLogical(0)));
1203              int lc = length(coef); double *cc = REAL(coef);              int lc = length(coef); double *cc = REAL(coef);
1204    
1205              ssclme_factor(x);              ssclme_factor(x);
1206              Rprintf("%3d %.3f", iter + 1, dev[REML ? 1 : 0]);      if (iter == 0) Rprintf("  EM iterations\n");
1207        Rprintf("%3d %.3f", iter, dev[REML ? 1 : 0]);
1208              for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);              for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1209              Rprintf("\n");              Rprintf("\n");
1210              UNPROTECT(1);              UNPROTECT(1);
1211          }          }
1212    
1213    /**
1214     * Perform ECME steps for the REML or ML criterion.
1215     *
1216     * @param x pointer to an ssclme object
1217     * @param nsteps pointer to an integer scalar - the number of ECME steps to perform
1218     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1219     * @param verb pointer to a logical scalar indicating verbose mode
1220     *
1221     * @return NULL
1222     */
1223    SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1224    {
1225        SEXP
1226            Omega = GET_SLOT(x, Matrix_OmegaSym);
1227        int
1228            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1229            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1230            REML = asLogical(REMLp),
1231            i, info, iter,
1232            nEM = asInteger(nsteps),
1233            nf = length(Omega),
1234            verbose = asLogical(verb);
1235        double
1236            *dev = REAL(GET_SLOT(x, Matrix_devianceSym));
1237    
1238        if (verbose) EMsteps_verbose_print(x, 0, REML);
1239        for (iter = 0; iter < nEM; iter++) {
1240            common_ECME_gradient(x, REML, Omega);
1241            status[0] = status[1] = 0;
1242            for (i = 0; i < nf; i++) {
1243                double *vali = REAL(VECTOR_ELT(Omega, i));
1244    
1245                F77_CALL(dpotrf)("U", nc + i, vali, nc + i, &info);
1246                if (info)
1247                    error("DPOTRF returned error code %d in Omega[[%d]] update",
1248                          info, i + 1);
1249                F77_CALL(dpotri)("U", nc + i, vali, nc + i, &info);
1250                if (info)
1251                    error("DPOTRI returned error code %d in Omega[[%d]] update",
1252                          info, i + 1);
1253            }
1254            if (verbose) EMsteps_verbose_print(x, iter + 1, REML);
1255      }      }
1256      ssclme_factor(x);      ssclme_factor(x);
1257      return R_NilValue;      return R_NilValue;
1258  }  }
1259    
1260  /**  /**
1261     * Evaluate the gradient with respect to the indicators of the
1262     * positions in the Omega matrices.
1263     *
1264     * @param x Pointer to an ssclme object
1265     * @param REML indicator of whether REML is to be used
1266     * @param val Pointer to an object of the same structure as Omega
1267     */
1268    static void indicator_gradient(SEXP x, int REML, SEXP val)
1269    {
1270        SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1271        int
1272            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1273            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1274            i, nf = length(Omega);
1275    
1276        common_ECME_gradient(x, REML, val);
1277        for (i = 0; i < nf; i++) {
1278            int info, nci = nc[i], mi = (Gp[i+1] - Gp[i])/nci;
1279            double
1280                *work = Memcpy((double *) Calloc(nci * nci, double),
1281                               REAL(VECTOR_ELT(Omega, i)), nci * nci),
1282                alpha = (double) -mi, beta = (double) mi;
1283    
1284            F77_CALL(dpotrf)("U", &nci, work, &nci, &info);
1285            if (info)
1286                error("DPOTRF gave error code %d on Omega[[%d]]", info, i + 1);
1287            F77_CALL(dtrtri)("U", "N", &nci, work, &nci, &info);
1288            if (info)
1289                error("DPOTRI gave error code %d on Omega[[%d]]", info, i + 1);
1290            F77_CALL(dsyrk)("U", "N", &nci, &nci, &alpha, work, &nci,
1291                            &beta, REAL(VECTOR_ELT(val, i)), &nci);
1292            Free(work);
1293        }
1294    }
1295    
1296    /**
1297     * Overwrite a gradient with respect to positions in Omega[[i]] by the
1298     * gradient with respect to the unconstrained parameters.
1299     *
1300     * @param grad pointer to a gradient wrt positions.  Contents are overwritten.
1301     * @param Omega pointer to a list of symmetric matrices (upper storage).
1302     */
1303    static void unconstrained_gradient(SEXP grad, SEXP Omega)
1304    {
1305        int i, ii, j, nf = length(Omega);
1306        double one = 1.0, zero = 0.0;
1307    
1308        for (i = 0; i < nf; i++) {
1309            SEXP Omegai = VECTOR_ELT(Omega, i);
1310            int nci = *INTEGER(getAttrib(Omegai, R_DimSymbol)),
1311                ncisq = nci * nci, ncip1 = nci + 1;
1312            double *chol =
1313                Memcpy(Calloc(ncisq, double), REAL(Omegai), ncisq),
1314                *ai = REAL(VECTOR_ELT(grad, i)),
1315                *tmp = Calloc(ncisq, double);
1316    
1317            F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1318            if (j)
1319                error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1320            /* calculate and store grad[i,,] %*% t(R) */
1321            for (j = 0; j < nci; j++) {
1322                F77_CALL(dsymv)("U", &nci, &one, ai, &nci, chol + j, &nci,
1323                                &zero, tmp + j, &nci);
1324            }
1325            /* full symmetric product gives diagonals */
1326            F77_CALL(dtrmm)("R", "U", "T", "N", &nci, &nci, &one, chol, &nci,
1327                            Memcpy(ai, tmp, ncisq), &nci);
1328            /* overwrite upper triangle with gradients for positions in L' */
1329            for (ii = 1; ii < nci; ii++) {
1330                for (j = 0; j < ii; j++) {
1331                    ai[j + ii*nci] = 2. * chol[j*ncip1] * tmp[j + ii*nci];
1332                    ai[ii + j*nci] = 0.;
1333                }
1334            }
1335            Free(chol); Free(tmp);
1336        }
1337    }
1338    
1339    /**
1340     * Fills cvec with unlist(lapply(mList, function(el) el[upper.tri(el, strict = FALSE)]))
1341     *
1342     * @param mList pointer to a list of REAL matrices
1343     * @param nc number of columns in each matrix
1344     * @param cvec pointer to REAL vector to receive the result
1345     */
1346    static void upperTriList_to_vector(SEXP mList, int *nc, SEXP cvec)
1347    {
1348        int i, ii, j, nf = length(mList), pos = 0;
1349    
1350        pos = 0;                    /* position in vector */
1351        for (i = 0; i < nf; i++) {
1352            SEXP omgi = VECTOR_ELT(mList, i);
1353            int nci = nc[i];
1354            for (j = 0; j < nci; j++) {
1355                for (ii = 0; ii <= j; ii++) {
1356                    REAL(cvec)[pos++] = REAL(omgi)[ii + j * nci];
1357                }
1358            }
1359        }
1360    }
1361    
1362    /**
1363   * Return the gradient of the ML or REML deviance.   * Return the gradient of the ML or REML deviance.
1364   *   *
1365   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1366   * @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
1367   * @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
1368     * @param Uncp pointer to a logical scalar indicating if result should be a single vector
1369   *   *
1370   * @return pointer to a numeric vector of the gradient.   * @return pointer to the gradient as a list of matrices or as a vector.
1371   */   */
1372    SEXP ssclme_grad(SEXP x, SEXP REMLp, SEXP Unc, SEXP OneVector)
1373    {
1374        SEXP ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1375    
1376        indicator_gradient(x, asLogical(REMLp), ans);
1377        if (asLogical(Unc))
1378            unconstrained_gradient(ans, GET_SLOT(x, Matrix_OmegaSym));
1379        if (asLogical(OneVector)) {
1380            int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));
1381            SEXP val = PROTECT(allocVector(REALSXP, coef_length(length(ans), nc)));
1382    
1383            upperTri_as_vector(ans, nc, val);
1384            UNPROTECT(2);
1385            return val;
1386        }
1387        UNPROTECT(1);
1388        return ans;
1389    }
1390    
1391  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1392  {  {
1393      SEXP      SEXP
# Line 1310  Line 1496 
1496  }  }
1497    
1498  /**  /**
1499     * Return the Hessian of the ML or REML deviance.  This is a
1500     * placeholder until I work out the evaluation of the analytic
1501     * Hessian, which probably will involve several helper functions.
1502     *
1503     * @param x pointer to an ssclme object
1504     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1505     * @param Uncp pointer to a logical scalar indicating if the
1506     * unconstrained parameterization is to be used
1507     *
1508     * @return pointer to an approximate Hessian matrix
1509     */
1510    SEXP ssclme_Hessian(SEXP x, SEXP REMLp, SEXP Uncp)
1511    {
1512        int j, ncoef = coef_length(length(GET_SLOT(x, Matrix_OmegaSym)),
1513                                   INTEGER(GET_SLOT(x, Matrix_ncSym))),
1514            unc = asLogical(Uncp);
1515        SEXP ans = PROTECT(allocMatrix(REALSXP, ncoef, ncoef)),
1516            base = PROTECT(ssclme_coef(x, Uncp)),
1517            current = PROTECT(duplicate(base)),
1518            gradient;
1519    
1520        for (j = 0; j < ncoef; j++) {
1521            double delta = (REAL(base)[j] ? 1.e-7 * REAL(base)[j] : 1.e-7);
1522            int i;
1523    
1524            for (i = 0; i < ncoef; i++) REAL(current)[i] = REAL(base)[i];
1525            REAL(current)[j] += delta/2.;
1526            ssclme_coefGets(x, current, Uncp);
1527            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1528            for (i = 0; i < ncoef; i++) REAL(ans)[j * ncoef + i] = REAL(gradient)[i];
1529            UNPROTECT(1);
1530            REAL(current)[j] -= delta;
1531            ssclme_coefGets(x, current, Uncp);
1532            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1533            for (i = 0; i < ncoef; i++)
1534                REAL(ans)[j * ncoef + i] =
1535                    (REAL(ans)[j * ncoef + i] - REAL(gradient)[i])/ delta;
1536            UNPROTECT(1);
1537            for (i = 0; i < j; i++) { /* symmetrize */
1538                REAL(ans)[j * ncoef + i] = REAL(ans)[i * ncoef + j] =
1539                    (REAL(ans)[j * ncoef + i] + REAL(ans)[i * ncoef + j])/2.;
1540            }
1541        }
1542        UNPROTECT(3);
1543        return ans;
1544    }
1545    
1546    /**
1547   * Calculate and return the fitted values.   * Calculate and return the fitted values.
1548   *   *
1549   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
# Line 1471  Line 1705 
1705      UNPROTECT(1);      UNPROTECT(1);
1706      return ans;      return ans;
1707  }  }
1708    

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

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