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 175, Fri May 28 00:28:19 2004 UTC revision 176, Fri May 28 13:23:44 2004 UTC
# Line 1  Line 1 
1  #include "ssclme.h"  #include "ssclme.h"
2    
3    #define slot_dup(dest, src, sym)  SET_SLOT(dest, sym, duplicate(GET_SLOT(src, sym)))
4    
5    /**
6     * Using the sscCrosstab object from the grouping factors, generate
7     * the slots in an ssclme object related to the symmetric sparse
8     * matrix representation of Z'Z.  If the model matrices for the
9     * grouping factors have only one column each then the structure can
10     * be copied, otherwise it must be generated from the sscCrosstab and
11     * the number of columns per grouping factor.
12     *
13     * @param nf number of factors
14     * @param nc vector of length nf+2 with number of columns in model matrices
15     * @param ctab pointer to the sscCrosstab object
16     * @param ssc pointer to an ssclme object to be filled out
17     */
18  static  static
19  void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)  void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)
20  {  {
21      int *snc, i, copyonly = 1;      int *snc, i, copyonly = 1;
22    
23      for (i = 0; i < nf; i++) {      SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));
24          if (nc[i] > 1) copyonly = 0;      snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
25        for (i = 0; i <= nf; i++) {
26            snc[i] = nc[i];
27            if (nc[i] > 1 && i < nf) copyonly = 0;
28      }      }
29      if (copyonly) {      if (copyonly) {
30          SET_SLOT(ssc, Matrix_pSym, duplicate(GET_SLOT(ctab, Matrix_pSym)));          slot_dup(ssc, ctab, Matrix_pSym);
31          SET_SLOT(ssc, Matrix_iSym, duplicate(GET_SLOT(ctab, Matrix_iSym)));          slot_dup(ssc, ctab, Matrix_iSym);
32          SET_SLOT(ssc, Matrix_xSym, duplicate(GET_SLOT(ctab, Matrix_xSym)));          slot_dup(ssc, ctab, Matrix_xSym);
33          SET_SLOT(ssc, Matrix_DimSym,          slot_dup(ssc, ctab, Matrix_DimSym);
34                   duplicate(GET_SLOT(ctab, Matrix_DimSym)));          slot_dup(ssc, ctab, Matrix_GpSym);
35          SET_SLOT(ssc, Matrix_GpSym, duplicate(GET_SLOT(ctab, Matrix_GpSym)));          return;
36      } else {      }
37        {
38          int          int
39              *GpIn = INTEGER(GET_SLOT(ctab, Matrix_GpSym)),              *GpIn = INTEGER(GET_SLOT(ctab, Matrix_GpSym)),
40              *GpOut,              *GpOut,
# Line 41  Line 60 
60              }              }
61          }          }
62          nOut = GpOut[nf];       /* size of output matrix */          nOut = GpOut[nf];       /* size of output matrix */
63          SET_SLOT(ssc, Matrix_DimSym,          SET_SLOT(ssc, Matrix_DimSym, allocVector(INTSXP, 2));
                  duplicate(GET_SLOT(ctab, Matrix_DimSym)));  
64          dims = INTEGER(GET_SLOT(ssc, Matrix_DimSym));          dims = INTEGER(GET_SLOT(ssc, Matrix_DimSym));
65          dims[0] = dims[1] = nOut;          dims[0] = dims[1] = nOut;
66          SET_SLOT(ssc, Matrix_pSym, allocVector(INTSXP, nOut + 1));          SET_SLOT(ssc, Matrix_pSym, allocVector(INTSXP, nOut + 1));
# Line 87  Line 105 
105                                                      * another diagonal element */                                                      * another diagonal element */
106                          int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];                          int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];
107                          Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);                          Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);
108                          AiOut[ApOut[mjj + 1] - 1] = mjj; /* maybe mjj-1? */                          AiOut[ApOut[mjj + 1] - 1] = mjj;
109                      }                      }
110                  }                  }
111              }              }
112          }          }
113          Free(map); Free(ncc);          Free(map); Free(ncc);
114      }      }
     SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));  
     snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));  
     for (i = 0; i <= nf; i++) {  
         snc[i] = nc[i];  
     }  
115  }  }
116    
117  /**  /**
# Line 120  Line 133 
133      Free(sz);      Free(sz);
134  }  }
135    
136    /**
137     * Create an ssclme object from a list of grouping factors, sorted in
138     * order of non-increasing numbers of levels, and an integer vector of
139     * the number of columns in the model matrices.  There is one more
140     * element in ncv than in facs.  The last element is the number of
141     * columns in the model matrix for the fixed effects plus the
142     * response.  (i.e. p+1)
143     *
144     * @param facs pointer to a list of grouping factors
145     * @param ncv pointer to an integer vector of number of columns per model matrix
146     *
147     * @return pointer to an ssclme object
148     */
149  SEXP  SEXP
150  ssclme_create(SEXP facs, SEXP ncv)  ssclme_create(SEXP facs, SEXP ncv)
151  {  {
# Line 216  Line 242 
242      return val;      return val;
243  }  }
244    
245    /**
246     * Copy information on Z'Z accumulated in the bVar array to Z'Z
247     *
248     * @param ncj number of columns in this block
249     * @param Gpj initial column for this group
250     * @param Gpjp initial column for the next group
251     * @param bVj pointer to the ncj x ncj x mj array to be filled
252     * @param Ap column pointer array for Z'Z
253     * @param Ai row indices for Z'Z
254     * @param Ax elements of Z'Z
255     */
256  static  static
257  void bVj_to_A(int ncj, int Gpj, int Gpjp, const double bVj[],  void bVj_to_A(int ncj, int Gpj, int Gpjp, const double bVj[],
258                const int Ap[], const int Ai[], double Ax[])                const int Ap[], const int Ai[], double Ax[])
# Line 232  Line 269 
269      }      }
270  }  }
271    
272    /**
273     * Copy the dimnames from the list of grouping factors and the model
274     * matrices for the grouping factors into the appropriate parts of the
275     * ssclme object.
276     *
277     * @param x pointer to an ssclme object
278     * @param facs pointer to a list of factors
279     * @param mmats pointer to a list of model matrices
280     *
281     * @return NULL
282     */
283  SEXP  SEXP
284  ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)  ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)
285  {  {
# Line 257  Line 305 
305      return R_NilValue;      return R_NilValue;
306  }  }
307    
308    /**
309     * Update the numerical entries x, ZtX, and XtX in an ssclme object
310     * according to a set of model matrices.
311     *
312     * @param x pointer to an ssclme object
313     * @param facs pointer to a list of grouping factors
314     * @param mmats pointer to a list of model matrices
315     *
316     * @return NULL
317     */
318  SEXP  SEXP
319  ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)  ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)
320  {  {
# Line 373  Line 431 
431      return R_NilValue;      return R_NilValue;
432  }  }
433    
434    /**
435     * Inflate Z'Z according to Omega and create the factorization LDL'
436     *
437     * @param x pointer to an ssclme object
438     *
439     * @return NULL
440     */
441  SEXP ssclme_inflate_and_factor(SEXP x)  SEXP ssclme_inflate_and_factor(SEXP x)
442  {  {
443      SEXP      SEXP
# Line 427  Line 492 
492      return R_NilValue;      return R_NilValue;
493  }  }
494    
495    
496    /**
497     * If status[["factored"]] is FALSE, create and factor Z'Z+Omega, then
498     * create RZX and RXX, the deviance components, and the value of the
499     * deviance for both ML and REML.
500     *
501     * @param x pointer to an ssclme object
502     *
503     * @return NULL
504     */
505  SEXP ssclme_factor(SEXP x)  SEXP ssclme_factor(SEXP x)
506  {  {
507      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
# Line 520  Line 595 
595      return R_NilValue;      return R_NilValue;
596  }  }
597    
598    /**
599     * Return the position of probe in the sorted index vector ind.  It is
600     * known that the position is greater than or equal to start so a linear
601     * search from start is used.
602     *
603     * @param probe value to be matched
604     * @param start index at which to start
605     * @param ind vector of indices
606     *
607     * @return index of the entry matching probe
608     */
609  static  static
610  int ldl_update_ind(int probe, int start, const int ind[])  int ldl_update_ind(int probe, int start, const int ind[])
611  {  {
# Line 529  Line 615 
615  }  }
616    
617  /**  /**
618   * Update the diagonal blocks of the inverse of LDL' (=Z'Z+W)   * Update the diagonal blocks of the inverse of LDL' (=Z'Z+W).  The
619     * lower Cholesky factors of the updated blocks are stored in the bVar
620     * slot.
621   *   *
622   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
623   *   *
# Line 620  Line 708 
708      return R_NilValue;      return R_NilValue;
709  }  }
710    
711    /**
712     * If necessary, factor Z'Z+Omega, ZtX, and XtX then, if necessary,
713     * form RZX, RXX, and bVar for the inverse of the Cholesky factor.
714     *
715     * @param x pointer to an ssclme object
716     *
717     * @return NULL (x is updated in place)
718     */
719  SEXP ssclme_invert(SEXP x)  SEXP ssclme_invert(SEXP x)
720  {  {
721      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
# Line 659  Line 755 
755      return R_NilValue;      return R_NilValue;
756  }  }
757    
758    /**
759     * Create and insert initial values for Omega_i.
760     *
761     * @param x pointer to an ssclme object
762     *
763     * @return NULL
764     */
765  SEXP ssclme_initial(SEXP x)  SEXP ssclme_initial(SEXP x)
766  {  {
767      SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),      SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),
# Line 783  Line 886 
886                                          nobs + 1 - pp1 : nobs))));                                          nobs + 1 - pp1 : nobs))));
887  }  }
888    
889    /**
890     * Calculate the length of the parameter vector, which is called coef
891     * for historical reasons.
892     *
893     * @param nf number of factors
894     * @param nc number of columns in the model matrices for each factor
895     *
896     * @return total length of the coefficient vector
897     */
898  static  static
899  int coef_length(int nf, const int nc[])  int coef_length(int nf, const int nc[])
900  {  {
# Line 792  Line 904 
904  }  }
905    
906  /**  /**
907   * Extract the upper triangles of the Omega matrices.   * Extract the upper triangles of the Omega matrices.  These aren't
908   * (These aren't "coefficients" but the extractor is   * "coefficients" but the extractor is called coef for historical
909   * called coef for historical reasons.)   * reasons.  Within each group these values are in the order of the
910     * diagonal entries first then the strict upper triangle in row
911     * order.
912   *   *
913   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
914   *   *
# Line 833  Line 947 
947    
948  /**  /**
949   * Extract the unconstrained parameters that determine the   * Extract the unconstrained parameters that determine the
950   * Omega matrices. (Called coef for historical reasons.)   * Omega matrices. (Called coef for historical reasons.)  The
951     * unconstrained parameters are derived from the LDL' decomposition of
952     * Omega_i.  The first nc[i] entries in each group are the diagonals
953     * of log(D) followed by the strict lower triangle of L in column
954     * order.
955   *   *
956   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
957   *   *
# Line 973  Line 1091 
1091      return x;      return x;
1092  }  }
1093    
1094    /**
1095     * Perform a number of ECME steps for the REML or ML criterion.
1096     *
1097     * @param x pointer to an ssclme object
1098     * @param nsteps pointer to an integer scalar giving the number of ECME steps to perform
1099     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1100     * @param verb pointer to a logical scalar indicating verbose mode
1101     *
1102     * @return NULL
1103     */
1104  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1105  {  {
1106      SEXP      SEXP
# Line 1065  Line 1193 
1193      return R_NilValue;      return R_NilValue;
1194  }  }
1195    
1196    /**
1197     * Return the gradient of the ML or REML deviance.
1198     *
1199     * @param x pointer to an ssclme object
1200     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1201     * @param Uncp pointer to a logical scalar indicating if the unconstrained parameterization is to be used
1202     *
1203     * @return pointer to a numeric vector of the gradient.
1204     */
1205  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1206  {  {
1207      SEXP      SEXP
1208          Omega = GET_SLOT(x, Matrix_OmegaSym),          Omega = GET_SLOT(x, Matrix_OmegaSym),
1209          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ncsl = GET_SLOT(x, Matrix_ncSym),  
1210          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1211      int      int
1212          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1213          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1214          *nc = INTEGER(ncsl),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1215          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1216          cind, i, n = dims[0],          cind, i, n = dims[0],
1217          nf = length(Omega),          nf = length(Omega),
# Line 1121  Line 1257 
1257          F77_CALL(dsyrk)("U", "N", &nci, &ki,          F77_CALL(dsyrk)("U", "N", &nci, &ki,
1258                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1259                          &alpha, tmp, &nci);                          &alpha, tmp, &nci);
1260          alpha = ((double)(REML?(nobs-p):nobs));          alpha = (double)(REML ? (nobs-p) : nobs);
1261          F77_CALL(dsyrk)("U", "N", &nci, &mi,          F77_CALL(dsyrk)("U", "N", &nci, &mi,
1262                          &alpha, b + Gp[i], &nci,                          &alpha, b + Gp[i], &nci,
1263                          &one, tmp, &nci);                          &one, tmp, &nci);
# Line 1139  Line 1275 
1275              int k, odind = cind + nci;              int k, odind = cind + nci;
1276              if (uncst) {              if (uncst) {
1277                  int ione = 1, kk;                  int ione = 1, kk;
1278                  double *rr = Calloc(nci, double);                  double *rr = Calloc(nci, double); /* j'th row of R, the Cholesky factor */
1279                  nlme_symmetrize(tmp, nci);                  nlme_symmetrize(tmp, nci);
1280                  for (j = 0; j < nci; j++, cind++) {                  for (j = 0; j < nci; j++, cind++) {
1281                      for (k = 0; k < nci; k++) rr[k] = chol[j + k*nci];                      for (k = 0; k < j; k++) rr[k] = 0.;
1282                        for (k = j; k < nci; k++) rr[k] = chol[j + k*nci];
1283                      REAL(ans)[cind] = 0.;                      REAL(ans)[cind] = 0.;
1284                      for (k = j; k < nci; k++) {                      for (k = j; k < nci; k++) {
1285                          for (kk = j; kk < nci; kk++) {                          for (kk = j; kk < nci; kk++) {
# Line 1150  Line 1287 
1287                                  tmp[kk * nci + k];                                  tmp[kk * nci + k];
1288                          }                          }
1289                      }                      }
                     for (k = 0; k < nci; k++) rr[k] *= rr[j];  
1290                      for (k = j + 1; k < nci; k++) {                      for (k = j + 1; k < nci; k++) {
1291                          REAL(ans)[odind++] =                          REAL(ans)[odind++] = 2. * rr[j] *
1292                              F77_CALL(ddot)(&nci, rr, &ione, tmp + k, &nci) +                              F77_CALL(ddot)(&nci, rr, &ione, tmp + k*nci, &ione);
                             F77_CALL(ddot)(&nci, rr, &ione,  
                                            tmp + k*nci, &ione);  
1293                      }                      }
1294                  }                  }
1295                  Free(rr);                  Free(rr);
# Line 1175  Line 1309 
1309      return ans;      return ans;
1310  }  }
1311    
1312    /**
1313     * Calculate and return the fitted values.
1314     *
1315     * @param x pointer to an ssclme object
1316     * @param facs list of grouping factors
1317     * @param mmats list of model matrices
1318     * @param useRf pointer to a logical scalar indicating if the random effects should be used
1319     *
1320     * @return pointer to a numeric array of fitted values
1321     */
1322  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats, SEXP useRf)  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats, SEXP useRf)
1323  {  {
1324      SEXP val, b;      SEXP val, b;
# Line 1251  Line 1395 
1395      return Omg;      return Omg;
1396  }  }
1397    
1398  #define slot_dup(sym)  SET_SLOT(ans, sym, duplicate(GET_SLOT(x, sym)))  /**
1399     * Copy an ssclme object collapsing the fixed effects slots to the response only.
1400     *
1401     * @param x pointer to an ssclme object
1402     *
1403     * @return a duplicate of x with the fixed effects slots collapsed to the response only
1404     */
1405  SEXP ssclme_collapse(SEXP x)  SEXP ssclme_collapse(SEXP x)
1406  {  {
1407      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),
# Line 1260  Line 1409 
1409          Dim = GET_SLOT(x, Matrix_DimSym);          Dim = GET_SLOT(x, Matrix_DimSym);
1410      int nf = length(Omega), nz = INTEGER(Dim)[1];      int nf = length(Omega), nz = INTEGER(Dim)[1];
1411    
1412      slot_dup(Matrix_DSym);      slot_dup(ans, x, Matrix_DSym);
1413      slot_dup(Matrix_DIsqrtSym);      slot_dup(ans, x, Matrix_DIsqrtSym);
1414      slot_dup(Matrix_DimSym);      slot_dup(ans, x, Matrix_DimSym);
1415      slot_dup(Matrix_GpSym);      slot_dup(ans, x, Matrix_GpSym);
1416      slot_dup(Matrix_LiSym);      slot_dup(ans, x, Matrix_LiSym);
1417      slot_dup(Matrix_LpSym);      slot_dup(ans, x, Matrix_LpSym);
1418      slot_dup(Matrix_LxSym);      slot_dup(ans, x, Matrix_LxSym);
1419      slot_dup(Matrix_OmegaSym);      slot_dup(ans, x, Matrix_OmegaSym);
1420      slot_dup(Matrix_ParentSym);      slot_dup(ans, x, Matrix_ParentSym);
1421      slot_dup(Matrix_bVarSym);      slot_dup(ans, x, Matrix_bVarSym);
1422      slot_dup(Matrix_devianceSym);      slot_dup(ans, x, Matrix_devianceSym);
1423      slot_dup(Matrix_devCompSym);      slot_dup(ans, x, Matrix_devCompSym);
1424      slot_dup(Matrix_iSym);      slot_dup(ans, x, Matrix_iSym);
1425      slot_dup(Matrix_ncSym);      slot_dup(ans, x, Matrix_ncSym);
1426      slot_dup(Matrix_statusSym);      slot_dup(ans, x, Matrix_statusSym);
1427      slot_dup(Matrix_pSym);      slot_dup(ans, x, Matrix_pSym);
1428      slot_dup(Matrix_xSym);      slot_dup(ans, x, Matrix_xSym);
1429      INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;      INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1430      SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));      SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));
1431      REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;      REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;
# Line 1289  Line 1438 
1438      return ans;      return ans;
1439  }  }
1440    
1441    
1442    /**
1443     * Create an lme object from its components.  This is not done by
1444     * new("lme", ...) at the R level because of the possibility of
1445     * causing the copying of very large objects.
1446     *
1447     * @param call Pointer to the original call
1448     * @param facs pointer to the list of grouping factors
1449     * @param x pointer to the model matrices (may be of length zero)
1450     * @param model pointer to the model frame
1451     * @param REML pointer to a logical scalar indicating if REML is used
1452     * @param rep pointer to the converged ssclme object
1453     * @param fitted pointer to the fitted values
1454     * @param residuals pointer to the residuals
1455     *
1456     * @return an lme object
1457     */
1458  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,
1459                     SEXP rep, SEXP fitted, SEXP residuals)                     SEXP rep, SEXP fitted, SEXP residuals)
1460  {  {

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

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