SCM

SCM Repository

[matrix] Diff of /branches/Matrix-mer2/src/lmer.c
ViewVC logotype

Diff of /branches/Matrix-mer2/src/lmer.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1091, Mon Dec 12 13:59:58 2005 UTC revision 1092, Tue Dec 13 00:26:11 2005 UTC
# Line 131  Line 131 
131    
132  #define flag_not_factored(x) *LOGICAL(GET_SLOT(x, Matrix_statusSym)) = 0  #define flag_not_factored(x) *LOGICAL(GET_SLOT(x, Matrix_statusSym)) = 0
133    
134    /**
135     * Create the diagonal blocks of the variance-covariance matrix of the
136     * random effects
137     *
138     * @param Linv - cholmod_sparse representation of L^{-1}
139     * @param nf - number of grouping factors
140     * @param Gp - group pointers
141     * @param nc - number of columns per factor
142     * @param bVar - list of 3-d arrays to be filled in
143     * @param uplo - "U" or "L" for upper or lower triangle
144     */
145  static void  static void
146  Linv_to_bVar(cholmod_sparse *Linv, int nf, const int Gp[], const int nc[],  Linv_to_bVar(cholmod_sparse *Linv, const int Gp[], const int nc[],
147               SEXP bVar, const char uplo[])               SEXP bVar, const char uplo[])
148  {  {
149      int *Lii = (int*)(Linv->i), *Lip = (int*)(Linv->p), i;      int *Lii = (int*)(Linv->i), *Lip = (int*)(Linv->p), i, nf = LENGTH(bVar);
150      double *Lix = (double*)(Linv->x), one[] = {1,0}, zero[] = {0,0};      double *Lix = (double*)(Linv->x), one[] = {1,0}, zero[] = {0,0};
151    
152      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
# Line 183  Line 194 
194  }  }
195    
196  /**  /**
197     * Evaluate the quadratic form in b defined by Omega^{-1}
198     *
199     * @param b vector of random effects
200     * @param Omega - list of dpoMatrix objects defining the pattern for Omega
201     * @param nf - number of grouping factors
202     * @param Gp - group pointers
203     * @param nc - number of columns per factor
204     *
205     * @return
206     */
207    static double
208    b_quadratic(const double b[], SEXP Omega, const int Gp[], const int nc[])
209    {
210        int i, ione = 1, nf = LENGTH(Omega);
211        double ans = 0., one[] = {1.,0.};
212    
213        for (i = 0; i < nf; i++) {
214            int nci = nc[i], ntot = Gp[i + 1] - Gp[i];
215            int nlev = ntot/nci;
216            double *bcp = Memcpy(Calloc(ntot, double), b + Gp[i], ntot),
217                *omgf = REAL(GET_SLOT(dpoMatrix_chol(VECTOR_ELT(Omega, i)), Matrix_xSym));
218    
219            F77_CALL(dtrsm)("L", "U", "T", "N", &nci, &nlev, one, omgf, &nci, bcp, &nci);
220            ans += F77_CALL(ddot)(&ntot, bcp, &ione, bcp, &ione);
221            Free(bcp);
222        }
223        return ans;
224    }
225    
226    /**
227   * Calculate the length of the parameter vector (historically called "coef"   * Calculate the length of the parameter vector (historically called "coef"
228   * even though these are not coefficients).   * even though these are not coefficients).
229   *   *
# Line 315  Line 356 
356      double tol;      /* convergence tolerance for IRLS iterations */      double tol;      /* convergence tolerance for IRLS iterations */
357  } glmer_struct, *GlmerStruct;  } glmer_struct, *GlmerStruct;
358    
359    static void
360    internal_mer_bVar(SEXP x)
361    {
362        int q = LENGTH(GET_SLOT(x, Matrix_rZySym));
363        cholmod_factor
364            *L = (cholmod_factor *) R_ExternalPtrAddr(GET_SLOT(x, Matrix_LSym));
365        cholmod_sparse *eye = cholmod_speye(q, q, CHOLMOD_REAL, &c), *Linv;
366        int *Perm = (int *)(L->Perm), *iperm = Calloc(q, int), i;
367    
368                                    /* Create Linv */
369        for (i = 0; i < q; i++) iperm[Perm[i]] = i; /* create the inverse permutation */
370                                    /* apply the inverse permutation to the identity matrix */
371        for (i = 0; i < q; i++) ((int*)(eye->i))[i] = iperm[i];
372        Linv = cholmod_spsolve(CHOLMOD_L, L, eye, &c);
373        cholmod_free_sparse(&eye, &c);
374        Linv_to_bVar(Linv, INTEGER(GET_SLOT(x, Matrix_GpSym)),
375                     INTEGER(GET_SLOT(x, Matrix_GpSym)),
376                     GET_SLOT(x, Matrix_bVarSym), "U");
377        cholmod_free_sparse(&Linv, &c);
378    }
379    
380  /**  /**
381   * Increment the fitted values from an mer object using the (possibly   * Increment the fitted values from an mer object using the (possibly
382   * externally stored) contents of the model matrices.   * externally stored) contents of the model matrices.
# Line 872  Line 934 
934      return ans;      return ans;
935  }  }
936    
937  static double Omega_log_det(SEXP Omega, int nf, int *nc, int *Gp)  static double
938    Omega_log_det(SEXP Omega, int nf, int *nc, int *Gp)
939  {  {
940      double ans = 0;      double ans = 0;
941      int i;      int i;
# Line 887  Line 950 
950      return ans;      return ans;
951  }  }
952    
953    /* FIXME: Separate the calculation of the offset random effects from
954     * the calculation of the deviance contributions. */
955    
956  /**  /**
957   * Determine the deviance components associated with each of the   * Determine the deviance components associated with each of the
958   * levels of a grouping factor at the conditional modes or a value   * levels of a grouping factor at the conditional modes or a value
959   * offset from the conditional modes by delb.   * offset from the conditional modes by delb.
960   *   *
961   * @param GS pointer to a GlmerStruct   * @param GS pointer to a GlmerStruct
962   * @param b pointer to the conditional modes of the random effects   * @param b conditional modes of the random effects
963   * @param nlev number of levels of the kth grouping factor   * @param Gp group pointers
964   * @param nc number of columns in the model matrix for the kth   * @param nc number of columns in the model matrix for the kth
965   * grouping factor   * grouping factor
966   * @param k index (0-based) of the grouping factor   * @param k index (0-based) of the grouping factor
# Line 904  Line 970 
970   * for this grouping factor   * for this grouping factor
971   * @param bVfac 3-dimensional array holding the factors of the   * @param bVfac 3-dimensional array holding the factors of the
972   * conditional variance-covariance matrix of the random effects   * conditional variance-covariance matrix of the random effects
973    FIXME: This is wrong.  It is bVar[[i]] not bVfac that is being passed.
974    This only affects the AGQ method.
975   * @param devcmp array to hold the deviance components   * @param devcmp array to hold the deviance components
976   *   *
977   * @return devcmp   * @return devcmp
978   */   */
979  static double*  static double*
980  rel_dev_1(GlmerStruct GS, SEXP b, int nlev, int nc, int k,  rel_dev_1(GlmerStruct GS, const double b[], int nlev, int nc, int k,
981            const double delb[], const double OmgFac[],            const double delb[], const double OmgFac[],
982            const double bVfac[], double devcmp[])            const double bVfac[], double devcmp[])
983  {  {
984      SEXP devs, flist = GET_SLOT(GS->mer, Matrix_flistSym);      SEXP devs;
985      int *fv = INTEGER(VECTOR_ELT(flist, k)), i, ione = 1,      int *fv = INTEGER(VECTOR_ELT(GET_SLOT(GS->mer, Matrix_flistSym), k)),
986          j, ntot = nlev * nc;          i, j;
987      double *bb = REAL(VECTOR_ELT(b, k)), *bcp = (double *) NULL;  /*     double *bcp = (double *) NULL; */
988    
989      AZERO(devcmp, nlev);      AZERO(devcmp, nlev);
990    #if 0
991      if (delb) {      if (delb) {
992            int ione = 1, ntot = nlev * nc;
993          double sumsq = 0;          double sumsq = 0;
994                                  /* copy the contents of b */                                  /* copy the contents of b */
995          bcp = Memcpy(Calloc(ntot, double), bb, ntot);          bcp = Memcpy(Calloc(ntot, double), b, ntot);
996          if (nc == 1) {          if (nc == 1) {
997              sumsq = delb[0] * delb[0];              sumsq = delb[0] * delb[0];
998              for (i = 0; i < nlev; i++) bb[i] += delb[0] * bVfac[i];              for (i = 0; i < nlev; i++) b[i] += delb[0] * bVfac[i];
999          } else {          } else {
1000              int ncsq = nc * nc;              int ncsq = nc * nc;
1001              double *tmp = Calloc(nc, double);              double *tmp = Calloc(nc, double);
# Line 933  Line 1003 
1003                  Memcpy(tmp, delb, nc);                  Memcpy(tmp, delb, nc);
1004                  F77_CALL(dtrmv)("U", "N", "N", &nc, &(bVfac[i * ncsq]),                  F77_CALL(dtrmv)("U", "N", "N", &nc, &(bVfac[i * ncsq]),
1005                                  &nc, tmp, &ione);                                  &nc, tmp, &ione);
1006                  for (j = 0; j < nc; j++) bb[i + j * nc] = tmp[j];                  for (j = 0; j < nc; j++) b[i + j * nc] = tmp[j];
1007              }              }
1008                                  /* sum of squares of delb */                                  /* sum of squares of delb */
1009              for (j = 0; j < nc; j++) sumsq += delb[j] * delb[j];              for (j = 0; j < nc; j++) sumsq += delb[j] * delb[j];
1010          }          }
1011          for (i = 0; i < nlev; i++) devcmp[i] = -sumsq;          for (i = 0; i < nlev; i++) devcmp[i] = -sumsq;
1012      }      }
1013    #endif
1014      Memcpy(REAL(GS->eta), GS->off, GS->n);      Memcpy(REAL(GS->eta), GS->off, GS->n);
1015  /*     fitted_ranef(flist, GS->unwtd, b, &nc, REAL(GS->eta)); */      internal_mer_fitted(
1016            GS->mer, (double *) NULL,
1017            REAL(GET_SLOT(GET_SLOT(GS->mer, Matrix_ZtSym), Matrix_xSym)),
1018            REAL(GS->eta));
1019      eval_check_store(GS->linkinv, GS->rho, GS->mu);      eval_check_store(GS->linkinv, GS->rho, GS->mu);
1020      devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));      devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));
1021      for (i = 0; i < GS->n; i++)      for (i = 0; i < GS->n; i++)
# Line 949  Line 1023 
1023      UNPROTECT(1);      UNPROTECT(1);
1024      if (nc == 1) {      if (nc == 1) {
1025          for (i = 0; i < nlev; i++) {          for (i = 0; i < nlev; i++) {
1026              double tmp = *OmgFac * bb[i];              double tmp = *OmgFac * b[i];
1027              devcmp[i] += tmp * tmp;              devcmp[i] += tmp * tmp;
1028          }          }
1029      } else {      } else {
# Line 957  Line 1031 
1031          int ione = 1;          int ione = 1;
1032    
1033          for (i = 0; i < nlev; i++) {          for (i = 0; i < nlev; i++) {
1034              for (j = 0; j < nc; j++) tmp[j] = bb[i + j * nlev];              for (j = 0; j < nc; j++) tmp[j] = b[i + j * nlev];
1035              F77_CALL(dtrmv)("U", "N", "N", &nc, OmgFac, &nc,              F77_CALL(dtrmv)("U", "N", "N", &nc, OmgFac, &nc,
1036                              tmp, &ione);                              tmp, &ione);
1037              for (j = 0; j < nc; j++)              for (j = 0; j < nc; j++)
1038                  devcmp[i] += tmp[j] * tmp[j];                  devcmp[i] += tmp[j] * tmp[j];
1039          }          }
1040      }      }
1041    #if 0
1042      if (delb) {      if (delb) {
1043          Memcpy(bb, bcp, ntot);          Memcpy(b, bcp, ntot);
1044          Free(bcp);          Free(bcp);
1045      }      }
1046    #endif
1047      return devcmp;      return devcmp;
1048  }  }
1049    
# Line 1021  Line 1097 
1097  }  }
1098    
1099  /**  /**
1100   * Evaluate the conditional deviance for a given set of fixed effects.   * Evaluate the conditional deviance for a given set of random effects.
1101   *   *
1102   * @param GS Pointer to a GlmerStruct   * @param GS Pointer to a GlmerStruct
1103   * @param fixed value of the fixed effects   * @param fixed value of the fixed effects
# Line 1029  Line 1105 
1105   * @return conditional deviance   * @return conditional deviance
1106   */   */
1107  static double  static double
1108  random_effects_deviance(GlmerStruct GS, SEXP b)  random_effects_deviance(GlmerStruct GS, const double b[])
1109  {  {
1110      SEXP devs;      SEXP devs;
1111      int i;      int i;
1112      double ans;      double ans;
1113    
1114      Memcpy(REAL(GS->eta), GS->off, GS->n);      Memcpy(REAL(GS->eta), GS->off, GS->n);
1115  /*     fitted_ranef(GET_SLOT(GS->mer, Matrix_flistSym), GS->unwtd, b, */      internal_mer_fitted(GS->mer, (double*) NULL, GS->Zt, REAL(GS->eta));
 /*               INTEGER(GET_SLOT(GS->mer, Matrix_ncSym)), REAL(GS->eta)); */  
1116      eval_check_store(GS->linkinv, GS->rho, GS->mu);      eval_check_store(GS->linkinv, GS->rho, GS->mu);
1117      devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));      devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));
1118      for (i = 0, ans = 0; i < GS->n; i++) ans += REAL(devs)[i];      for (i = 0, ans = 0; i < GS->n; i++) ans += REAL(devs)[i];
# Line 1045  Line 1120 
1120      return ans;      return ans;
1121  }  }
1122    
1123    /**
1124     * Return the components of the deviance corresponding to the levels
1125     * of a single grouping factor.
1126     *
1127     * @param GS
1128     * @param f
1129     *
1130     * @return
1131     */
1132    static double*
1133    glmer_dev_comp_1d(GlmerStruct GS, double f[])
1134    {
1135        SEXP devs;
1136        int i, *fl1 = INTEGER(VECTOR_ELT(GET_SLOT(GS->mer, Matrix_flistSym), 0));
1137    
1138        Memcpy(REAL(GS->eta), GS->off, GS->n);
1139        eval_check_store(GS->linkinv, GS->rho, GS->mu);
1140        devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));
1141        for (i = 0; i < GS->n; i++) f[fl1[i] - 1] += REAL(devs)[i];
1142        return f;
1143    }
1144    
1145    
1146  static void  static void
1147  internal_glmer_ranef_update(GlmerStruct GS, SEXP b)  internal_glmer_ranef_update(GlmerStruct GS, SEXP b)
1148  {  {
# Line 1052  Line 1150 
1150          bVar = GET_SLOT(GS->mer, Matrix_bVarSym),          bVar = GET_SLOT(GS->mer, Matrix_bVarSym),
1151          Omega = GET_SLOT(GS->mer, Matrix_OmegaSym);          Omega = GET_SLOT(GS->mer, Matrix_OmegaSym);
1152      int i, ione = 1, j, k;      int i, ione = 1, j, k;
1153    /*     double *b = REAL(GET_SLOT(GS->mer, Matrix_ranefSym); */
1154      double devr, one = 1;      double devr, one = 1;
1155    
1156  /*     bhat = PROTECT(lmer_ranef(GS->mer)); */      bhat = R_NilValue;
1157                                  /* subtract deviance at b */                                  /* subtract deviance at b */
1158      devr = -random_effects_deviance(GS, b);  /*     devr = -random_effects_deviance(GS, b); */
1159      for (i = 0; i < GS->nf; i++) {      for (i = 0; i < GS->nf; i++) {
1160          SEXP Bi = VECTOR_ELT(b, i);          SEXP Bi = VECTOR_ELT(b, i);
1161          int *dims = INTEGER(getAttrib(Bi, R_DimSymbol));          int *dims = INTEGER(getAttrib(Bi, R_DimSymbol));
# Line 1073  Line 1172 
1172                               REAL(VECTOR_ELT(Omega, i)),                               REAL(VECTOR_ELT(Omega, i)),
1173                               ncisqr);                               ncisqr);
1174    
1175                                  /* subtract quadratic form in                                  /* subtract quadratic form in */
1176                                   * Omega[[i]] at b  */                                  /* Omega[[i]] at b  */
1177          F77_CALL(dpotrf)("U", &nci, omgfac, &nci, &j);          F77_CALL(dpotrf)("U", &nci, omgfac, &nci, &j);
1178          if (j)          if (j)
1179              error(_("Leading %d minor of Omega[[%d]] not positive definite"),              error(_("Leading %d minor of Omega[[%d]] not positive definite"),
# Line 1107  Line 1206 
1206              for (j = 0; j < nci; j++)              for (j = 0; j < nci; j++)
1207                  bpropi[k + j * nlev] = bhati[k + j * nlev] + delta[j];                  bpropi[k + j * nlev] = bhati[k + j * nlev] + delta[j];
1208          }          }
1209                                  /* add quadratic form in                                  /* add quadratic form in */
1210                                   * Omega[[i]] at bprop  */                                  /* Omega[[i]] at bprop  */
1211          Memcpy(bcopy, bpropi, ntot);          Memcpy(bcopy, bpropi, ntot);
1212          F77_CALL(dtrmm)("R", "U", "T", "N", &nlev, &nci, &one,          F77_CALL(dtrmm)("R", "U", "T", "N", &nlev, &nci, &one,
1213                          omgfac, &nci, bcopy, &nlev);                          omgfac, &nci, bcopy, &nlev);
# Line 1117  Line 1216 
1216          Free(bcopy); Free(chol); Free(delta); Free(omgfac);          Free(bcopy); Free(chol); Free(delta); Free(omgfac);
1217      }      }
1218                                  /* add deviance at bprop */                                  /* add deviance at bprop */
1219      devr += random_effects_deviance(GS, bprop);  /*     devr += random_effects_deviance(GS, bprop); */
1220    
1221      if (unif_rand() < exp(-0.5 * devr))      if (unif_rand() < exp(-0.5 * devr))
1222          for (i = 0; i < GS->nf; i++) { /* copy each face of b */          for (i = 0; i < GS->nf; i++) { /* copy each face of b */
# Line 1394  Line 1493 
1493   *   *
1494   * @return the approximation to the deviance as computed using AGQ   * @return the approximation to the deviance as computed using AGQ
1495   */   */
1496    SEXP glmer_devLaplace(SEXP pars, SEXP GSp)
1497    {
1498        GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);
1499        SEXP bVar = GET_SLOT(GS->mer, Matrix_bVarSym),
1500            Omega = GET_SLOT(GS->mer, Matrix_OmegaSym);
1501        int *Gp = INTEGER(GET_SLOT(GS->mer, Matrix_GpSym)),
1502            *nc = INTEGER(GET_SLOT(GS->mer, Matrix_ncSym)),
1503            i, j, k;
1504        double *bhat = REAL(GET_SLOT(GS->mer, Matrix_ranefSym));
1505    
1506        if (!isReal(pars) || LENGTH(pars) != GS->npar)
1507            error(_("`%s' must be a numeric vector of length %d"),
1508                  "pars", GS->npar);
1509        if (!internal_bhat(GS, REAL(pars), REAL(pars) + (GS->p)))
1510            return ScalarReal(DBL_MAX);
1511        internal_mer_bVar(GS->mer);
1512        return ScalarReal(Omega_log_det(Omega, LENGTH(Omega), nc, Gp) +
1513                          random_effects_deviance(GS, bhat) +
1514                          b_quadratic(bhat, Omega, Gp, nc));
1515    }
1516    
1517    /**
1518     * Compute the approximation to the deviance using adaptive
1519     * Gauss-Hermite quadrature (AGQ).  When nAGQ == 1 this is the Laplace
1520     * approximation.
1521     *
1522     * @param pars pointer to a numeric vector of parameters
1523     * @param GSp pointer to a GlmerStruct object
1524     * @param nAGQp pointer to a scalar integer representing the number of
1525     * points in AGQ to use
1526     *
1527     * @return the approximation to the deviance as computed using AGQ
1528     */
1529  SEXP glmer_devAGQ(SEXP pars, SEXP GSp, SEXP nAGQp)  SEXP glmer_devAGQ(SEXP pars, SEXP GSp, SEXP nAGQp)
1530  {  {
1531      GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);      GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);
1532      SEXP Omega, bVar;      SEXP Omega = GET_SLOT(GS->mer, Matrix_OmegaSym),
1533            bVar = GET_SLOT(GS->mer, Matrix_bVarSym);
1534      int i, j, k, nAGQ = asInteger(nAGQp);      int i, j, k, nAGQ = asInteger(nAGQp);
1535      int n2 = (nAGQ + 1)/2;      int n2 = (nAGQ + 1)/2,
1536            *Gp = INTEGER(GET_SLOT(GS->mer, Matrix_GpSym)),
1537            *nc = INTEGER(GET_SLOT(GS->mer, Matrix_ncSym));
1538      double *f0, LaplaceDev = 0, AGQadjst = 0,      double *f0, LaplaceDev = 0, AGQadjst = 0,
1539          *bhat = REAL(GET_SLOT(GS->mer, Matrix_ranefSym));          *bhat = REAL(GET_SLOT(GS->mer, Matrix_ranefSym));
1540    
# Line 1412  Line 1547 
1547      }      }
1548      if (!internal_bhat(GS, REAL(pars), REAL(pars) + (GS->p)))      if (!internal_bhat(GS, REAL(pars), REAL(pars) + (GS->p)))
1549          return ScalarReal(DBL_MAX);          return ScalarReal(DBL_MAX);
     return R_NilValue;  
     bVar = GET_SLOT(GS->mer, Matrix_bVarSym);  
     Omega = GET_SLOT(GS->mer, Matrix_OmegaSym);  
1550    
1551      for (i = 0; i < GS->nf; i++) {      for (i = 0; i < GS->nf; i++) {
1552          int *dims = INTEGER(getAttrib(VECTOR_ELT(bVar, i),          int nci = nc[i];
1553                                        R_DimSymbol));          int ncip1 = nci + 1, ncisqr = nci * nci,
1554          int nci = dims[0], nlev = dims[2];              nlev = (Gp[i + 1] - Gp[i])/nci;
1555          int ncip1 = nci + 1, ncisqr = nci * nci;          double *omgf = REAL(GET_SLOT(dpoMatrix_chol(VECTOR_ELT(Omega, i)), Matrix_xSym)),
1556          double *omg = Memcpy(Calloc(ncisqr, double),              *bVi = Memcpy(Calloc(ncisqr * nlev, double),
                              REAL(VECTOR_ELT(Omega, i)), ncisqr),  
             *bvar = Memcpy(Calloc(ncisqr * nlev, double),  
1557                             REAL(VECTOR_ELT(bVar, i)), ncisqr * nlev);                             REAL(VECTOR_ELT(bVar, i)), ncisqr * nlev);
1558    
         LaplaceDev = 0;  
                                 /* Calculate difference of determinants */  
         F77_CALL(dpotrf)("U", &nci, omg, &nci, &j);  
         if (j)  
             error(_("Leading %d minor of Omega[[%d]] not positive definite"),  
                   j, i + 1);  
1559          for (j = 0; j < nci; j++) { /* nlev * logDet(Omega_i) */          for (j = 0; j < nci; j++) { /* nlev * logDet(Omega_i) */
1560              /* Note: we subtract because Omega has been inverted */              LaplaceDev += 2 * nlev * log(omgf[j * ncip1]);
             LaplaceDev -= 2 * nlev * log(omg[j * ncip1]);  
1561          }          }
1562          for (k = 0; k < nlev; k++) {          for (k = 0; k < nlev; k++) {
1563              double *bVik = bvar + k * ncisqr;              double *bVik = bVi + k * ncisqr;
1564              F77_CALL(dpotrf)("U", &nci, bVik, &nci, &j);              F77_CALL(dpotrf)("U", &nci, bVik, &nci, &j);
1565              if (j)              if (j)
1566                  error(_("Leading %d minor of bVar[[%d]][,,%d] not positive definite"),                  error(_("Leading %d minor of bVar[[%d]][,,%d] not positive definite"),
# Line 1447  Line 1570 
1570    
1571          f0 = Calloc(nlev, double);          f0 = Calloc(nlev, double);
1572          rel_dev_1(GS, bhat, nlev, nci, i, (double *) NULL,          rel_dev_1(GS, bhat, nlev, nci, i, (double *) NULL,
1573                    omg, bvar, f0);                    omgf, bVi, f0);
1574          for (k = 0; k < nlev; k++) LaplaceDev += f0[k];          for (k = 0; k < nlev; k++) LaplaceDev += f0[k];
1575          if (nAGQ > 1) {          if (nAGQ > 1) {
1576              double *fx = Calloc(nlev, double),              double *fx = Calloc(nlev, double),
# Line 1460  Line 1583 
1583                  delb[0] = GHQ_x[nAGQ][k];                  delb[0] = GHQ_x[nAGQ][k];
1584                  if (delb[0]) {                  if (delb[0]) {
1585                      rel_dev_1(GS, bhat, nlev, nci, i, delb,                      rel_dev_1(GS, bhat, nlev, nci, i, delb,
1586                                omg, bvar, fx);                                omgf, bVi, fx);
1587                      for (j = 0; j < nlev; j++) {                      for (j = 0; j < nlev; j++) {
1588                          rellik[j] += GHQ_w[nAGQ][k] *                          rellik[j] += GHQ_w[nAGQ][k] *
1589                              exp(-(fx[j] - f0[j])/2);                              exp(-(fx[j] - f0[j])/2);
1590                      }                      }
1591                      delb[0] *= -1;                      delb[0] *= -1;
1592                      rel_dev_1(GS, bhat, nlev, nci, i, delb,                      rel_dev_1(GS, bhat, nlev, nci, i, delb,
1593                                omg, bvar, fx);                                omgf, bVi, fx);
1594                      for (j = 0; j < nlev; j++) {                      for (j = 0; j < nlev; j++) {
1595                          rellik[j] += GHQ_w[nAGQ][k] *                          rellik[j] += GHQ_w[nAGQ][k] *
1596                              exp(-(fx[j] - f0[j])/2);                              exp(-(fx[j] - f0[j])/2);
# Line 1481  Line 1604 
1604                  AGQadjst -= 2 * log(rellik[j]);                  AGQadjst -= 2 * log(rellik[j]);
1605              Free(fx); Free(rellik);              Free(fx); Free(rellik);
1606          }          }
1607          Free(f0); Free(omg); Free(bvar);          Free(f0); Free(bVi);
1608      }      }
     UNPROTECT(1);  
1609      return ScalarReal(LaplaceDev + AGQadjst);      return ScalarReal(LaplaceDev + AGQadjst);
1610  }  }
1611    
# Line 2329  Line 2451 
2451          int q = LENGTH(ranefP), p = LENGTH(GET_SLOT(x, Matrix_rXySym));          int q = LENGTH(ranefP), p = LENGTH(GET_SLOT(x, Matrix_rXySym));
2452          cholmod_factor          cholmod_factor
2453              *L = (cholmod_factor *) R_ExternalPtrAddr(GET_SLOT(x, Matrix_LSym));              *L = (cholmod_factor *) R_ExternalPtrAddr(GET_SLOT(x, Matrix_LSym));
         cholmod_sparse *eye = cholmod_speye(q, q, CHOLMOD_REAL, &c), *Linv;  
2454          cholmod_dense *RZX = as_cholmod_dense(RZXP), *tmp1;          cholmod_dense *RZX = as_cholmod_dense(RZXP), *tmp1;
2455          int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
2456              *Perm = (int *)(L->Perm),              *Perm = (int *)(L->Perm),
# Line 2355  Line 2476 
2476          F77_CALL(dtrsm)("R", "U", "N", "N", &q, &p, m1,          F77_CALL(dtrsm)("R", "U", "N", "N", &q, &p, m1,
2477                          REAL(GET_SLOT(GET_SLOT(x, Matrix_RXXSym), Matrix_xSym)),                          REAL(GET_SLOT(GET_SLOT(x, Matrix_RXXSym), Matrix_xSym)),
2478                          &p, RZXinv, &q);                          &p, RZXinv, &q);
                                 /* Create Linv */  
         /* apply the inverse permutation to the identity matrix */  
         for (i = 0; i < q; i++) ((int *)(eye->i))[i] = iperm[i];  
         Linv = cholmod_spsolve(CHOLMOD_L, L, eye, &c);  
         cholmod_free_sparse(&eye, &c);  
         Linv_to_bVar(Linv, nf, Gp, nc, bVarP, "U");  
2479    
2480            internal_mer_bVar(x);
2481          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
2482              int nci = nc[i], RZXrows = Gp[i + 1] - Gp[i];              int nci = nc[i], RZXrows = Gp[i + 1] - Gp[i];
2483              int ncisq = nci * nci, nlev = RZXrows/nci;              int ncisq = nci * nci, nlev = RZXrows/nci;
# Line 2409  Line 2525 
2525              }              }
2526              Free(tmp);              Free(tmp);
2527          }          }
         cholmod_free_sparse(&Linv, &c);  
2528          Free(iperm);          Free(iperm);
2529          status[2] = 1; status[3] = 0;          status[2] = 1; status[3] = 0;
2530      }      }

Legend:
Removed from v.1091  
changed lines
  Added in v.1092

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