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 245, Mon Jul 12 12:46:54 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 1091  Line 1091 
1091      return x;      return x;
1092  }  }
1093    
1094    
1095    /**
1096     * Returns the inverse of the updated Omega matrices for an ECME
1097     * iteration.  These matrices are also used in the gradient calculation.
1098     *
1099     * @param x pointer to an ssclme object
1100     * @param REML indicator of REML being used
1101     * @param val pointer to a list of symmetric q_i by q_i matrices
1102     */
1103    static void
1104    common_ECME_gradient(SEXP x, int REML, SEXP val)
1105    {
1106        SEXP
1107            RZXsl = GET_SLOT(x, Matrix_RZXSym),
1108            bVar = GET_SLOT(x, Matrix_bVarSym);
1109        int
1110            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1111            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1112            i, j, n = *INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1113            nf = length(val), nobs = nc[nf + 1], p = nc[nf] - 1;
1114        double
1115            *RZX = REAL(RZXsl),
1116            *b = RZX + p * INTEGER(getAttrib(RZXsl, R_DimSymbol))[0],
1117            one = 1.0, zero = 0.0;
1118    
1119        ssclme_invert(x);
1120        for (i = 0; i < nf; i++) {
1121            int ki = Gp[i+1] - Gp[i],
1122                nci = nc[i],
1123                mi = ki/nci;
1124            double
1125                *vali = REAL(VECTOR_ELT(val, i)),
1126                alpha = (double)(REML ? (nobs-p) : nobs)/((double) mi);
1127    
1128            F77_CALL(dsyrk)("U", "N", &nci, &mi,
1129                            &alpha, b + Gp[i], &nci,
1130                            &zero, vali, &nci);
1131            alpha = 1./((double) mi);
1132            F77_CALL(dsyrk)("U", "N", &nci, &ki,
1133                            &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1134                            &one, vali, &nci);
1135            if (REML) {
1136                for (j = 0; j < p; j++) {
1137                    F77_CALL(dsyrk)("U", "N", &nci, &mi,
1138                                    &alpha, RZX + Gp[i] + j*n, &nci,
1139                                    &one, vali, &nci);
1140                }
1141            }
1142        }
1143    }
1144    
1145  /**  /**
1146   * Perform a number of ECME steps for the REML or ML criterion.   * Perform a number of ECME steps for the REML or ML criterion.
1147   *   *
# Line 1104  Line 1155 
1155  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1156  {  {
1157      SEXP      SEXP
1158          Omega = GET_SLOT(x, Matrix_OmegaSym),          Omega = GET_SLOT(x, Matrix_OmegaSym);
         RZXsl = GET_SLOT(x, Matrix_RZXSym),  
         ncsl = GET_SLOT(x, Matrix_ncSym),  
         bVar = GET_SLOT(x, Matrix_bVarSym);  
1159      int      int
1160          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
         *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),  
         *nc = INTEGER(ncsl),  
1161          *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),          *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1162          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1163          i, info, iter,          i, info, iter,
         n = dims[0],  
1164          nEM = asInteger(nsteps),          nEM = asInteger(nsteps),
1165          nf = length(ncsl) - 2,          nf = length(Omega),
         nobs = nc[nf + 1],  
         p,  
         pp1 = dims[1],  
1166          verbose = asLogical(verb);          verbose = asLogical(verb);
1167      double      double
1168          *RZX = REAL(RZXsl),          *dev = REAL(GET_SLOT(x, Matrix_devianceSym));
         *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),  
         *b,  
         alpha,  
         one = 1.,  
         zero = 0.;  
1169    
     p = pp1 - 1;  
     b = RZX + p * n;  
1170      if (verbose) {      if (verbose) {
1171          SEXP coef = PROTECT(ssclme_coef(x));          SEXP coef = PROTECT(ssclme_coef(x));
1172          int lc = length(coef); double *cc = REAL(coef);          int lc = length(coef); double *cc = REAL(coef);
# Line 1144  Line 1179 
1179          UNPROTECT(1);          UNPROTECT(1);
1180      }      }
1181      for (iter = 0; iter < nEM; iter++) {      for (iter = 0; iter < nEM; iter++) {
1182          ssclme_invert(x);          common_ECME_gradient(x, REML, Omega);
1183            status[0] = status[1] = 0;
1184          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
1185              int ki = Gp[i+1] - Gp[i],              double *vali = REAL(VECTOR_ELT(Omega, i));
                 nci = nc[i],  
                 mi = ki/nci;  
             double  
                 *vali = REAL(VECTOR_ELT(Omega, i));  
1186    
1187              alpha = ((double)(REML?(nobs-p):nobs))/((double)mi);              F77_CALL(dpotrf)("U", nc + i, vali, nc + i, &info);
             F77_CALL(dsyrk)("U", "N", &nci, &mi,  
                             &alpha, b + Gp[i], &nci,  
                             &zero, vali, &nci);  
             alpha = 1./((double) mi);  
             F77_CALL(dsyrk)("U", "N", &nci, &ki,  
                             &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,  
                             &one, vali, &nci);  
             if (REML) {  
                 int j;  
                 for (j = 0; j < p; j++) {  
                     F77_CALL(dsyrk)("U", "N", &nci, &mi,  
                                 &alpha, RZX + Gp[i] + j*n, &nci,  
                                 &one, vali, &nci);  
                 }  
             }  
             F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);  
1188              if (info)              if (info)
1189                  error("DPOTRF returned error code %d in Omega[[%d]] update",                  error("DPOTRF returned error code %d in Omega[[%d]] update",
1190                        info, i + 1);                        info, i + 1);
1191              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);              F77_CALL(dpotri)("U", nc + i, vali, nc + i, &info);
1192              if (info)              if (info)
1193                  error("DPOTRI returned error code %d in Omega[[%d]] update",                  error("DPOTRI returned error code %d in Omega[[%d]] update",
1194                        info, i + 1);                        info, i + 1);
1195          }          }
         status[0] = status[1] = 0;  
1196          if (verbose) {          if (verbose) {
1197              SEXP coef = PROTECT(ssclme_coef(x));              SEXP coef = PROTECT(ssclme_coef(x));
1198              int lc = length(coef); double *cc = REAL(coef);              int lc = length(coef); double *cc = REAL(coef);
# Line 1194  Line 1209 
1209  }  }
1210    
1211  /**  /**
1212     * Evaluate the gradient with respect to the indicators of the
1213     * positions in the Omega matrices.
1214     *
1215     * @param x Pointer to an ssclme object
1216     * @param REML indicator of whether REML is to be used
1217     * @param val Pointer to an object of the same structure as Omega
1218     */
1219    static void indicator_gradient(SEXP x, int REML, SEXP val)
1220    {
1221        SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1222        int
1223            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1224            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1225            i, nf = length(Omega);
1226    
1227        common_ECME_gradient(x, REML, val);
1228        for (i = 0; i < nf; i++) {
1229            int info, nci = nc[i], mi = (Gp[i+1] - Gp[i])/nci;
1230            double
1231                *work = Memcpy((double *) Calloc(nci * nci, double),
1232                               REAL(VECTOR_ELT(Omega, i)), nci * nci),
1233                alpha = (double) -mi, beta = (double) mi;
1234    
1235            F77_CALL(dpotrf)("U", &nci, work, &nci, &info);
1236            if (info)
1237                error("DPOTRF gave error code %d on Omega[[%d]]", info, i + 1);
1238            F77_CALL(dtrtri)("U", "N", &nci, work, &nci, &info);
1239            if (info)
1240                error("DPOTRI gave error code %d on Omega[[%d]]", info, i + 1);
1241            F77_CALL(dsyrk)("U", "N", &nci, &nci, &alpha, work, &nci,
1242                            &beta, REAL(VECTOR_ELT(val, i)), &nci);
1243            Free(work);
1244        }
1245    }
1246    
1247    /**
1248     * Overwrite a gradient with respect to positions in Omega[[i]] by the
1249     * gradient with respect to the unconstrained parameters.
1250     *
1251     * @param grad pointer to a gradient wrt positions.  Overwritten.
1252     * @param Omega pointer to a list of symmetric (upper storage)
1253     * matrices Omega[[i]].
1254     */
1255    static void unconstrained_gradient(SEXP grad, SEXP Omega)
1256    {
1257        int i, ii, j, nf = length(Omega);
1258        double one = 1.0, zero = 0.0;
1259    
1260        for (i = 0; i < nf; i++) {
1261            SEXP Omegai = VECTOR_ELT(Omega, i);
1262            int nci = *INTEGER(getAttrib(Omegai, R_DimSymbol)),
1263                ncisq = nci * nci, ncip1 = nci + 1;
1264            double *chol =
1265                Memcpy(Calloc(ncisq, double), REAL(Omegai), ncisq),
1266                *ai = REAL(VECTOR_ELT(grad, i)),
1267                *tmp = Calloc(ncisq, double);
1268    
1269            F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1270            if (j)
1271                error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1272            /* calculate and store grad[i,,] %*% t(R) */
1273            for (j = 0; j < nci; j++) {
1274                F77_CALL(dsymv)("U", &nci, &one, ai, &nci, chol + j, &nci,
1275                                &zero, tmp + j, &nci);
1276            }
1277            /* full symmetric product gives diagonals */
1278            F77_CALL(dtrmm)("R", "U", "T", "N", &nci, &nci, &one, chol, &nci,
1279                            Memcpy(ai, tmp, ncisq), &nci);
1280            /* overwrite lower triangle with gradients for positions in L; */
1281            for (ii = 1; ii < nci; ii++) {
1282                for (j = 0; j < ii; j++) {
1283                    ai[ii + j*nci] = 2. * chol[j*ncip1] * tmp[j + ii*nci];
1284                    ai[j + ii*nci] = 0.;
1285                }
1286            }
1287            Free(chol); Free(tmp);
1288        }
1289    }
1290    
1291    SEXP ssclme_grad(SEXP x, SEXP REMLp, SEXP Uncst)
1292    {
1293        SEXP ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1294    
1295        indicator_gradient(x, asLogical(REMLp), ans);
1296        if (asLogical(Uncst))
1297            unconstrained_gradient(ans, GET_SLOT(x, Matrix_OmegaSym));
1298        UNPROTECT(1);
1299        return ans;
1300    }
1301    
1302    
1303    /**
1304   * Return the gradient of the ML or REML deviance.   * Return the gradient of the ML or REML deviance.
1305   *   *
1306   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
# Line 1310  Line 1417 
1417  }  }
1418    
1419  /**  /**
1420     * Return the Hessian of the ML or REML deviance.  This is a
1421     * placeholder until I work out the evaluation of the analytic
1422     * Hessian, which probably will involve several helper functions.
1423     *
1424     * @param x pointer to an ssclme object
1425     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1426     * @param Uncp pointer to a logical scalar indicating if the
1427     * unconstrained parameterization is to be used
1428     *
1429     * @return pointer to an approximate Hessian matrix
1430     */
1431    SEXP ssclme_Hessian(SEXP x, SEXP REMLp, SEXP Uncp)
1432    {
1433        int j, ncoef = coef_length(length(GET_SLOT(x, Matrix_OmegaSym)),
1434                                   INTEGER(GET_SLOT(x, Matrix_ncSym))),
1435            unc = asLogical(Uncp);
1436        SEXP ans = PROTECT(allocMatrix(REALSXP, ncoef, ncoef)),
1437            base = PROTECT(unc ? ssclme_coefUnc(x) : ssclme_coef(x)),
1438            current = PROTECT(duplicate(base)),
1439            gradient;
1440    
1441        for (j = 0; j < ncoef; j++) {
1442            double delta = (REAL(base)[j] ? 1.e-7 * REAL(base)[j] : 1.e-7);
1443            int i;
1444    
1445            for (i = 0; i < ncoef; i++) REAL(current)[i] = REAL(base)[i];
1446            REAL(current)[j] += delta/2.;
1447            if (unc) {
1448                ssclme_coefGetsUnc(x, current);
1449            } else {
1450                ssclme_coefGets(x, current);
1451            }
1452            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1453            for (i = 0; i < ncoef; i++) REAL(ans)[j * ncoef + i] = REAL(gradient)[i];
1454            UNPROTECT(1);
1455            REAL(current)[j] -= delta;
1456            if (unc) {
1457                ssclme_coefGetsUnc(x, current);
1458            } else {
1459                ssclme_coefGets(x, current);
1460            }
1461            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1462            for (i = 0; i < ncoef; i++)
1463                REAL(ans)[j * ncoef + i] = (REAL(ans)[j * ncoef + i] - REAL(gradient)[i])/
1464                    delta;
1465            UNPROTECT(1);
1466            /* symmetrize */
1467            for (i = 0; i < j; i++) {
1468                REAL(ans)[j * ncoef + i] = REAL(ans)[i * ncoef + j] =
1469                    (REAL(ans)[j * ncoef + i] + REAL(ans)[i * ncoef + j])/2.;
1470            }
1471        }
1472        UNPROTECT(3);
1473        return ans;
1474    }
1475    
1476    /**
1477   * Calculate and return the fitted values.   * Calculate and return the fitted values.
1478   *   *
1479   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
# Line 1471  Line 1635 
1635      UNPROTECT(1);      UNPROTECT(1);
1636      return ans;      return ans;
1637  }  }
1638    

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

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