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 28, Sun Mar 28 17:06:05 2004 UTC revision 125, Sun Apr 25 11:59:37 2004 UTC
# Line 52  Line 52 
52   *   *
53   * @return pointer to an integer R vector.   * @return pointer to an integer R vector.
54   */   */
55  static  
56  SEXP ctab_permute(SEXP ctab)  SEXP ctab_permute(SEXP ctab)
57  {  {
58      SEXP val, GpSl = GET_SLOT(ctab, Matrix_GpSym);      SEXP val, GpSl = GET_SLOT(ctab, Matrix_GpSym);
# Line 65  Line 65 
65          j,          j,
66          n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],          n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
67          nf = length(GpSl) - 1,          nf = length(GpSl) - 1,
         nz = Ap[n],             /* number of non-zeros */  
68          pos;          pos;
69    
70      if (ctab_isNested(n, nf, 1, Ap, Ai, Gp))      if (ctab_isNested(n, nf, 1, Ap, Ai, Gp))
# Line 73  Line 72 
72      val =  allocVector(INTSXP, n);      val =  allocVector(INTSXP, n);
73      perm = INTEGER(val);      perm = INTEGER(val);
74      work = (int *) R_alloc(n, sizeof(int));      work = (int *) R_alloc(n, sizeof(int));
75      ssc_metis_order(n, nz, Ap, Ai, work, perm); /* perm gets inverse perm */      ssc_metis_order(n, Ap, Ai, work, perm);     /* perm gets inverse perm */
76      /* work now contains desired permutation but with groups scrambled */      /* work now contains desired permutation but with groups scrambled */
77    
78      /* copy work into perm preserving the order of the groups */      /* copy work into perm preserving the order of the groups */
# Line 189  Line 188 
188      }      }
189  }  }
190    
 static  
191  void ssclme_fill_LIp(int n, const int Parent[], int LIp[])  void ssclme_fill_LIp(int n, const int Parent[], int LIp[])
192  {  {
193      int *sz = Calloc(n, int), i;      int *sz = Calloc(n, int), i;
# Line 254  Line 252 
252      SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));      SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));
253      SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));      SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));
254      SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));      SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));
255          /* Zero the symmetric matrices (for cosmetic reasons only). */                                  /* Zero symmetric matrices (cosmetic) */
256      memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,      memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,
257             sizeof(double) * pp1 * pp1);             sizeof(double) * pp1 * pp1);
258      memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,      memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,
# Line 434  Line 432 
432          for (k = j+1; k < nf; k++) { /* off-diagonals */          for (k = j+1; k < nf; k++) { /* off-diagonals */
433              int *fpk = INTEGER(VECTOR_ELT(facs, k)),              int *fpk = INTEGER(VECTOR_ELT(facs, k)),
434                  *Apk = Ap + Gp[k],                  *Apk = Ap + Gp[k],
435                  nck = nc[k];                  nck = nc[k],
436                    scalar = ncj == 1 && nck == 1;
437              double              double
438                  *Zk = Z[k];                  *Zk = Z[k], *work;
439                if (!scalar) work = Calloc(ncj * nck, double);
440              for (i = 0; i < nobs; i++) {              for (i = 0; i < nobs; i++) {
441                  int ii, ind = -1, fpji = fpj[i] - 1,                  int ii, ind = -1, fpji = fpj[i] - 1,
442                      row = Gp[j] + fpji * ncj,                      row = Gp[j] + fpji * ncj,
443                      fpki = fpk[i] - 1,                      fpki = fpk[i] - 1,
444                      lastind = Apk[fpki + 1];                      lastind = Apk[fpki*nck + 1];
445                  for (ii = Apk[fpki]; ii < lastind; ii++) {                  for (ii = Apk[fpki*nck]; ii < lastind; ii++) {
446                      if (Ai[ii] == row) {                      if (Ai[ii] == row) {
447                          ind = ii;                          ind = ii;
448                          break;                          break;
449                      }                      }
450                  }                  }
451                  if (ind < 0) error("logic error in ssclme_update_mm");                  if (ind < 0) error("logic error in ssclme_update_mm");
452                  if (Ncj || nck > 1) {                  if (scalar) {   /* update scalars directly */
453                                  /* FIXME: run a loop to update */                      Ax[ind] += Zj[i] * Zk[i];
454                      error("code not yet written");                  } else {
455                  } else {        /* update scalars directly */                      int jj, offset = ind - Apk[fpki * nck];
456                      Ax[ind] += Zj[fpji] * Zk[fpki];                      F77_CALL(dgemm)("T", "N", &ncj, &nck, &ione, &one,
457                                        Zj + i, &nobs, Zk + i, &nobs,
458                                        &zero, work, &ncj);
459                        for (jj = 0; jj < nck; jj++) {
460                            ind = Apk[fpki * nck + jj] + offset;
461                            if (Ai[ind] != row)
462                                error("logic error in ssclme_update_mm");
463                            for (ii = 0; ii < ncj; ii++) {
464                                Ax[ind++] += work[jj * ncj + ii];
465                            }
466                  }                  }
467              }              }
468          }          }
469                if (!scalar) Free(work);
470            }
471      }      }
472      Free(Z);      Free(Z);
473      ssclme_transfer_dimnames(x, facs, mmats);      ssclme_transfer_dimnames(x, facs, mmats);
474      return R_NilValue;      return R_NilValue;
475  }  }
476    
477  SEXP ssclme_inflate_and_factor(SEXP lme)  SEXP ssclme_inflate_and_factor(SEXP x)
478  {  {
479      SEXP      SEXP
480          GpSlot = GET_SLOT(lme, Matrix_GpSym),          GpSlot = GET_SLOT(x, Matrix_GpSym),
481          Omega = GET_SLOT(lme, Matrix_OmegaSym);          Omega = GET_SLOT(x, Matrix_OmegaSym);
482      int n = INTEGER(GET_SLOT(lme, Matrix_DimSym))[1];      int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
483      int      int
484          *Ai = INTEGER(GET_SLOT(lme, Matrix_iSym)),          *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
485          *Ap = INTEGER(GET_SLOT(lme, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
486          *Flag = Calloc(n, int),          *Flag = Calloc(n, int),
487          *Gp = INTEGER(GpSlot),          *Gp = INTEGER(GpSlot),
488          *Lnz = Calloc(n, int),          *Lnz = Calloc(n, int),
489          *Pattern = Calloc(n, int),          *Pattern = Calloc(n, int),
490          *nc = INTEGER(GET_SLOT(lme, Matrix_ncSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
491          j,          j,
492          nf = length(GpSlot) - 1;          nf = length(GpSlot) - 1;
493      double      double
494          *D = REAL(GET_SLOT(lme, Matrix_DSym)),          *D = REAL(GET_SLOT(x, Matrix_DSym)),
495          *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),          *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
496          *Y = Calloc(n, double),          *Y = Calloc(n, double),
497          *xcp = Calloc(Ap[n], double);          *xcp = Calloc(Ap[n], double);
498    
499      Memcpy(xcp, REAL(GET_SLOT(lme, Matrix_xSym)), Ap[n]);      Memcpy(xcp, REAL(GET_SLOT(x, Matrix_xSym)), Ap[n]);
500      for (j = 0; j < nf; j++) {      for (j = 0; j < nf; j++) {
501          int  diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];          int  diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];
502          double *omgj = REAL(VECTOR_ELT(Omega, j));          double *omgj = REAL(VECTOR_ELT(Omega, j));
# Line 504  Line 514 
514          }          }
515      }      }
516      j = ldl_numeric(n, Ap, Ai, xcp,      j = ldl_numeric(n, Ap, Ai, xcp,
517                      INTEGER(GET_SLOT(lme, Matrix_LpSym)),                      INTEGER(GET_SLOT(x, Matrix_LpSym)),
518                      INTEGER(GET_SLOT(lme, Matrix_ParentSym)),                      INTEGER(GET_SLOT(x, Matrix_ParentSym)),
519                      Lnz, INTEGER(GET_SLOT(lme, Matrix_LiSym)),                      Lnz, INTEGER(GET_SLOT(x, Matrix_LiSym)),
520                      REAL(GET_SLOT(lme, Matrix_LxSym)),                      REAL(GET_SLOT(x, Matrix_LxSym)),
521                      D, Y, Pattern, Flag,                      D, Y, Pattern, Flag,
522                      (int *) NULL, (int *) NULL); /* P & Pinv */                      (int *) NULL, (int *) NULL); /* P & Pinv */
523      if (j != n)      if (j != n)
# Line 518  Line 528 
528      return R_NilValue;      return R_NilValue;
529  }  }
530    
531  SEXP ssclme_factor(SEXP lme)  SEXP ssclme_factor(SEXP x)
532  {  {
533      int *status = LOGICAL(GET_SLOT(lme, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
534    
535      if (!status[0]) {      if (!status[0]) {
536          SEXP          SEXP
537              GpSlot = GET_SLOT(lme, Matrix_GpSym),              GpSlot = GET_SLOT(x, Matrix_GpSym),
538              Omega = GET_SLOT(lme, Matrix_OmegaSym);              Omega = GET_SLOT(x, Matrix_OmegaSym);
539          int          int
540              *Gp = INTEGER(GpSlot),              *Gp = INTEGER(GpSlot),
541              *Li = INTEGER(GET_SLOT(lme, Matrix_LiSym)),              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
542              *Lp = INTEGER(GET_SLOT(lme, Matrix_LpSym)),              *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
543              *nc = INTEGER(GET_SLOT(lme, Matrix_ncSym)),              *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
544              i,              i,
545              n = INTEGER(GET_SLOT(lme, Matrix_DimSym))[1],              n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1],
546              nf = length(GpSlot) - 1,              nf = length(GpSlot) - 1,
547              nobs = nc[nf + 1],              nobs = nc[nf + 1],
548              nreml = nobs + 1 - nc[nf],              nreml = nobs + 1 - nc[nf],
549              pp1 = nc[nf],              pp1 = nc[nf],
550              pp2 = pp1 + 1;              pp2 = pp1 + 1;
551          double          double
552              *D = REAL(GET_SLOT(lme, Matrix_DSym)),              *D = REAL(GET_SLOT(x, Matrix_DSym)),
553              *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),              *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
554              *Lx = REAL(GET_SLOT(lme, Matrix_LxSym)),              *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
555              *RXX = REAL(GET_SLOT(lme, Matrix_RXXSym)),              *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
556              *RZX = REAL(GET_SLOT(lme, Matrix_RZXSym)),              *RZX = REAL(GET_SLOT(x, Matrix_RZXSym)),
557              *dcmp = REAL(getAttrib(lme, Matrix_devCompSym)),              *dcmp = REAL(getAttrib(x, Matrix_devCompSym)),
558              *deviance = REAL(getAttrib(lme, Matrix_devianceSym)),              *deviance = REAL(getAttrib(x, Matrix_devianceSym)),
559              minus1 = -1.,              minus1 = -1.,
560              one = 1.;              one = 1.;
561    
562          ssclme_inflate_and_factor(lme);          ssclme_inflate_and_factor(x);
563                                  /* Accumulate logdet of ZtZ+W */                                  /* Accumulate logdet of ZtZ+W */
564          dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;          dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;
565          for (i = 0; i < n; i++) dcmp[0] += log(D[i]);          for (i = 0; i < n; i++) dcmp[0] += log(D[i]);
# Line 579  Line 589 
589              }              }
590          }          }
591                                  /* ldl_lsolve on Z'X */                                  /* ldl_lsolve on Z'X */
592          Memcpy(RZX, REAL(GET_SLOT(lme, Matrix_ZtXSym)), n * pp1);          Memcpy(RZX, REAL(GET_SLOT(x, Matrix_ZtXSym)), n * pp1);
593          for (i = 0; i < pp1; i++) {          for (i = 0; i < pp1; i++) {
594              int j;              int j;
595              double *RZXi = RZX + i * n;              double *RZXi = RZX + i * n;
# Line 587  Line 597 
597              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
598          }          }
599                                  /* downdate and factor X'X */                                  /* downdate and factor X'X */
600          Memcpy(RXX, REAL(GET_SLOT(lme, Matrix_XtXSym)), pp1 * pp1);          Memcpy(RXX, REAL(GET_SLOT(x, Matrix_XtXSym)), pp1 * pp1);
601          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
602                          RZX, &n, &one, RXX, &pp1);                          RZX, &n, &one, RXX, &pp1);
603          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
604          if (i)          if (i) {
605              error("DPOTRF returned error code %d", i);              warning("Could not factor downdated X'X, code %d", i);
606                dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;
607            } else {
608                                  /* logdet of RXX */                                  /* logdet of RXX */
609          for (i = 0; i < (pp1 - 1); i++)          for (i = 0; i < (pp1 - 1); i++)
610              dcmp[2] += 2 * log(RXX[i*pp2]);              dcmp[2] += 2 * log(RXX[i*pp2]);
# Line 602  Line 614 
614              dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));              dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));
615          deviance[1] = dcmp[0] - dcmp[1] + /* REML */          deviance[1] = dcmp[0] - dcmp[1] + /* REML */
616              dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));              dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));
617            }
618          status[0] = 1;          status[0] = 1;
619          status[1] = 0;          status[1] = 0;
620      }      }
# Line 625  Line 638 
638   * @return R_NilValue (x is updated in place)   * @return R_NilValue (x is updated in place)
639    
640   */   */
641    static
642  SEXP ldl_inverse(SEXP x)  SEXP ldl_inverse(SEXP x)
643  {  {
644      SEXP      SEXP
# Line 697  Line 711 
711              }              }
712              tmp = Calloc(nr * nci, double); /* scratch storage */              tmp = Calloc(nr * nci, double); /* scratch storage */
713              nr1 = nr + 1;              nr1 = nr + 1;
714                                  /* initialize bVi to zero (cosmetic) */                                  /* initialize bVi to zero */
715              memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);              memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
716              for (j = G1; j < G2; j += nci) {              for (j = G1; j < G2; j += nci) {
717                  memset(tmp, 0, sizeof(double) * nr * nci);                  memset(tmp, 0, sizeof(double) * nr * nci);
# Line 713  Line 727 
727                          tmp[k * nr + kk] *= DIsqrt[LIi[LIp[j] + kk - 1]];                          tmp[k * nr + kk] *= DIsqrt[LIi[LIp[j] + kk - 1]];
728                      }                      }
729                  }                  }
730                  F77_CALL(dsyrk)("U", "T", &nci, &rr, &one, tmp, &nr,                  F77_CALL(dsyrk)("L", "T", &nci, &rr, &one, tmp, &nr,
731                                  &zero, bVi + (j - G1) * nci, &nci);                                  &zero, bVi + (j - G1) * nci, &nci);
732                  F77_CALL(dpotrf)("U", &nci, bVi + (j - G1) * nci,                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1) * nci,
733                                   &nci, &kk);                                   &nci, &kk);
734                  if (kk)         /* should never happen */                  if (kk)         /* should never happen */
735                      error(                      error(
736                          "Rank deficient variance matrix at group %d, level %d",                          "Rank deficient variance matrix at group %d, level %d",
737                          i + 1, j + 1);                          i + 1, j + 1);
738              }              }
739                Free(tmp);
740          }          }
741          return R_NilValue;          return R_NilValue;
742      }      }
# Line 778  Line 793 
793                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];
794                      }                      }
795                  }                  }
796                  F77_CALL(dsyrk)("U", "T", &nci, &nr, &one, tmp, &nr,                  F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,
797                                  &zero, mpt + (j - Gp[i])*nci, &nci);                                  &zero, mpt + (j - Gp[i])*nci, &nci);
798                  F77_CALL(dpotrf)("U", &nci, mpt + (j - Gp[i])*nci,                  F77_CALL(dpotrf)("L", &nci, mpt + (j - Gp[i])*nci,
799                                   &nci, &info);                                   &nci, &info);
800                  if (info)       /* should never happen */                  if (info)       /* should never happen */
801                      error(                      error(
# Line 793  Line 808 
808      return R_NilValue;      return R_NilValue;
809  }  }
810    
811  SEXP ssclme_invert(SEXP lme)  SEXP ssclme_invert(SEXP x)
812  {  {
813      int *status = LOGICAL(GET_SLOT(lme, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
814      if (!status[0]) ssclme_factor(lme);      if (!status[0]) ssclme_factor(x);
815        if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))
816            error("Unable to invert singular factor of downdated X'X");
817      if (!status[1]) {      if (!status[1]) {
818          SEXP          SEXP
819              RZXsl = GET_SLOT(lme, Matrix_RZXSym);              RZXsl = GET_SLOT(x, Matrix_RZXSym);
820          int          int
821              *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),              *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
822              *Li = INTEGER(GET_SLOT(lme, Matrix_LiSym)),              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
823              *Lp = INTEGER(GET_SLOT(lme, Matrix_LpSym)),              *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
824              i,              i,
825              n = dims[0],              n = dims[0],
826              pp1 = dims[1];              pp1 = dims[1];
827          double          double
828              *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),              *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
829              *Lx = REAL(GET_SLOT(lme, Matrix_LxSym)),              *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
830              *RXX = REAL(GET_SLOT(lme, Matrix_RXXSym)),              *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
831              *RZX = REAL(RZXsl),              *RZX = REAL(RZXsl),
832              one = 1.;              one = 1.;
833    
# Line 824  Line 841 
841              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
842              ldl_ltsolve(n, RZXi, Lp, Li, Lx);              ldl_ltsolve(n, RZXi, Lp, Li, Lx);
843          }          }
844          ldl_inverse(lme);          ldl_inverse(x);
845          status[1] = 1;          status[1] = 1;
846      }      }
847      return R_NilValue;      return R_NilValue;
# Line 869  Line 886 
886    
887  /**  /**
888   * Extract the conditional estimates of the fixed effects   * Extract the conditional estimates of the fixed effects
  * FIXME: Add names  
889   *   *
890   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
891   *   *
# Line 897  Line 913 
913    
914  /**  /**
915   * Extract the conditional modes of the random effects.   * Extract the conditional modes of the random effects.
  * FIXME: Change the returned value to be a named list of matrices  
  *        with dimnames.  
916   *   *
917   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
918   *   *
# Line 924  Line 938 
938      ssclme_invert(x);      ssclme_invert(x);
939      ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];      ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];
940      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
941          int nci = nc[i], Mi = (Gp[i+1] - Gp[i]), mi = Mi/nci;          int nci = nc[i], Mi = Gp[i+1] - Gp[i];
942          double *mm;          double *mm;
943    
944          SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, mi));          SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, Mi/nci));
945          mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);          mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);
946          b += nci * mi;          b += Mi;
947          for (j = 0; j < Mi; j++) mm[j] /= ryyinv;          for (j = 0; j < Mi; j++) mm[j] /= ryyinv;
948      }      }
949      UNPROTECT(1);      UNPROTECT(1);
950      return val;      return val;
951  }  }
952    
953  /**  /**
954   * Extract the ML or REML conditional estimate of sigma   * Extract the ML or REML conditional estimate of sigma
955   *   *
# Line 947  Line 962 
962  {  {
963      SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);      SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
964      int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],      int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],
965          nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))          nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))[
966          [length(GET_SLOT(x, Matrix_GpSym))];              length(GET_SLOT(x, Matrix_GpSym))];
967    
968      ssclme_invert(x);      ssclme_invert(x);
969      return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *      return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *
# Line 966  Line 981 
981    
982  /**  /**
983   * Extract the upper triangles of the Omega matrices.   * Extract the upper triangles of the Omega matrices.
984   * (These are not in any sense "coefficients" but the extractor is   * (These aren't "coefficients" but the extractor is
985   * called coef for historical reasons.)   * called coef for historical reasons.)
986   *   *
987   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
# Line 984  Line 999 
999    
1000      vind = 0;      vind = 0;
1001      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1002          int j, k, nci = nc[i];          int nci = nc[i];
1003            if (nci == 1) {
1004                vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];
1005            } else {
1006                int j, k, odind = vind + nci, ncip1 = nci + 1;
1007          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
1008    
1009          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
1010              for (k = 0; k <= j; k++) {                  vv[vind++] = omgi[j * ncip1];
1011                  vv[vind++] = omgi[j*nci + k];                  for (k = j + 1; k < nci; k++) {
1012                        vv[odind++] = omgi[k*nci + j];
1013              }              }
1014          }          }
1015                vind = odind;
1016            }
1017      }      }
1018      UNPROTECT(1);      UNPROTECT(1);
1019      return val;      return val;
1020  }  }
1021    
1022  /**  /**
1023   * Extract the upper triangles of the Omega matrices in the unconstrained   * Extract the unconstrained parameters that determine the
1024   * parameterization.   * Omega matrices. (Called coef for historical reasons.)
  * (These are not in any sense "coefficients" but the extractor is  
  * called coef for historical reasons.)  
1025   *   *
1026   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1027   *   *
1028   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of unconstrained parameters that determine the
1029   * Omega matrices   * Omega matrices
1030   */   */
1031  SEXP ssclme_coefUnc(SEXP x)  SEXP ssclme_coefUnc(SEXP x)
# Line 1026  Line 1047 
1047                                   REAL(VECTOR_ELT(Omega, i)), ncisq);                                   REAL(VECTOR_ELT(Omega, i)), ncisq);
1048              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1049              if (j)              /* should never happen */              if (j)              /* should never happen */
1050                  error("DPOTRF returned error code %d", j);                  error("DPOTRF returned error code %d on Omega[[%d]]",
1051                          j, i+1);
1052              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1053                  double diagj = tmp[j * ncip1];                  double diagj = tmp[j * ncip1];
1054                  vv[vind++] = 2. * log(diagj);                  vv[vind++] = 2. * log(diagj);
# Line 1047  Line 1069 
1069  }  }
1070    
1071  /**  /**
1072   * Assign the upper triangles of the Omega matrices in the unconstrained   * Assign the Omega matrices from the unconstrained parameterization.
  * parameterization.  
1073   *   *
1074   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1075   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1079  Line 1100 
1100              double              double
1101                  *omgi = REAL(VECTOR_ELT(Omega, i)),                  *omgi = REAL(VECTOR_ELT(Omega, i)),
1102                  *tmp = Calloc(ncisq, double),                  *tmp = Calloc(ncisq, double),
1103                  diagj, one = 1.;                  diagj, one = 1., zero = 0.;
1104              /* FIXEME: Change this to use a factor and dsyrk */  
                                 /* LD in omgi and L' in tmp */  
1105              memset(omgi, 0, sizeof(double) * ncisq);              memset(omgi, 0, sizeof(double) * ncisq);
1106              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1107                  omgi[j * ncip1] = diagj = exp(cc[cind++]);                  tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1108                  for (k = j + 1; k < nci; k++) {                  for (k = j + 1; k < nci; k++) {
1109                      omgi[j*nci + k] = diagj * (tmp[k*nci + j] = cc[odind++]);                      tmp[k*nci + j] = cc[odind++] * diagj;
1110                  }                  }
1111              }              }
1112              F77_CALL(dtrmm)("R", "U", "N", "U", &nci, &nci, &one,              F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1113                              tmp, &nci, omgi, &nci);                              tmp, &nci, &zero, omgi, &nci);
1114              Free(tmp);              Free(tmp);
1115              cind = odind;              cind = odind;
1116          }          }
# Line 1101  Line 1121 
1121    
1122  /**  /**
1123   * Assign the upper triangles of the Omega matrices.   * Assign the upper triangles of the Omega matrices.
1124   * (These are not in any sense "coefficients" but are   * (Called coef for historical reasons.)
  * called coef for historical reasons.)  
1125   *   *
1126   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1127   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1122  Line 1141 
1141                coef_length(nf, nc));                coef_length(nf, nc));
1142      cind = 0;      cind = 0;
1143      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1144          int j, k, nci = nc[i];          int nci = nc[i];
1145            if (nci == 1) {
1146                REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];
1147            } else {
1148                int j, k, odind = cind + nci, ncip1 = nci + 1;
1149          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
1150    
1151          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
1152              for (k = 0; k <= j; k++) {                  omgi[j * ncip1] = cc[cind++];
1153                  omgi[j*nci + k] = cc[cind++];                  for (k = j + 1; k < nci; k++) {
1154                        omgi[k*nci + j] = cc[odind++];
1155              }              }
1156          }          }
1157                cind = odind;
1158            }
1159      }      }
1160      status[0] = status[1] = 0;      status[0] = status[1] = 0;
1161      return x;      return x;
# Line 1192  Line 1219 
1219                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1220                              &one, vali, &nci);                              &one, vali, &nci);
1221              if (REML) {              if (REML) {
1222                  int mp = mi * p;                  int j;
1223                  F77_CALL(dsyrk)("U", "N", &nci, &mp,                  for (j = 0; j < p; j++) {
1224                                  &alpha, RZX + Gp[i], &nci,                      F77_CALL(dsyrk)("U", "N", &nci, &mi,
1225                                    &alpha, RZX + Gp[i] + j*n, &nci,
1226                                  &one, vali, &nci);                                  &one, vali, &nci);
1227              }              }
1228                }
1229              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1230              if (info) error("DPOTRF returned error code %d", info);              if (info)
1231                    error("DPOTRF returned error code %d in Omega[[%d]] update",
1232                          info, i + 1);
1233              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1234              if (info) error("DPOTRF returned error code %d", info);              if (info)
1235                    error("DPOTRI returned error code %d in Omega[[%d]] update",
1236                          info, i + 1);
1237          }          }
1238          status[0] = status[1] = 0;          status[0] = status[1] = 0;
1239      }      }
# Line 1208  Line 1241 
1241      return R_NilValue;      return R_NilValue;
1242  }  }
1243    
1244  SEXP ssclme_asSscMatrix(SEXP x)  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1245  {  {
1246      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix")));      SEXP
1247      int *dims = INTEGER(GET_SLOT(val, Matrix_DimSym));          Omega = GET_SLOT(x, Matrix_OmegaSym),
1248            RZXsl = GET_SLOT(x, Matrix_RZXSym),
1249            ncsl = GET_SLOT(x, Matrix_ncSym),
1250            bVar = GET_SLOT(x, Matrix_bVarSym);
1251        int
1252            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1253            *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1254            *nc = INTEGER(ncsl),
1255            REML = asLogical(REMLp),
1256            cind, i, n = dims[0],
1257            nf = length(Omega),
1258            nobs, p, pp1 = dims[1],
1259            uncst = asLogical(Uncp);
1260        double
1261            *RZX = REAL(RZXsl),
1262            *b,
1263            alpha,
1264            one = 1.;
1265        SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1266        int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1267    
1268      dims[0] = dims[1] = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];      ssclme_factor(x);
1269      SET_SLOT(val, Matrix_pSym, duplicate(GET_SLOT(x, Matrix_pSym)));      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {
1270      SET_SLOT(val, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_iSym)));          int ncoef = coef_length(nf, nc);
1271      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));          for (i = 0; i < ncoef; i++) REAL(ans)[i] = NA_REAL;
     CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0))[0] = 'U';  
1272      UNPROTECT(1);      UNPROTECT(1);
1273      return val;          return ans;
1274  }  }
1275        nobs = nc[nf + 1];
1276        p = pp1 - 1;
1277        b = RZX + p * n;
1278        ssclme_invert(x);
1279        cind = 0;
1280        for (i = 0; i < nf; i++) {
1281            int j, ki = Gp[i+1] - Gp[i],
1282                nci = nc[i], ncip1 = nci + 1, ncisq = nci * nci,
1283                mi = ki/nci;
1284            double
1285                *chol = Memcpy(Calloc(ncisq, double),
1286                               REAL(VECTOR_ELT(Omega, i)), ncisq),
1287                *tmp = Calloc(ncisq, double);
1288    
 SEXP ssclme_asTscMatrix(SEXP x)  
 {  
     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix")));  
     int *dims = INTEGER(GET_SLOT(val, Matrix_DimSym));  
1289    
1290      dims[0] = dims[1] = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1291      SET_SLOT(val, Matrix_pSym, duplicate(GET_SLOT(x, Matrix_LpSym)));          if (j)
1292      SET_SLOT(val, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_LiSym)));              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1293      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_LxSym)));          Memcpy(tmp, chol, ncisq);
1294      CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0))[0] = 'U';          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1295      CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0] = 'U';          if (j)
1296                error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
1297            alpha = (double) -mi;
1298            F77_CALL(dsyrk)("U", "N", &nci, &ki,
1299                            &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1300                            &alpha, tmp, &nci);
1301            alpha = ((double)(REML?(nobs-p):nobs));
1302            F77_CALL(dsyrk)("U", "N", &nci, &mi,
1303                            &alpha, b + Gp[i], &nci,
1304                            &one, tmp, &nci);
1305            if (REML) {
1306                for (j = 0; j < p; j++) {
1307                    F77_CALL(dsyrk)("U", "N", &nci, &mi,
1308                                    &one, RZX + Gp[i] + j*n, &nci,
1309                                    &one, tmp, &nci);
1310                }
1311            }
1312            if (nci == 1) {
1313                REAL(ans)[cind++] = *tmp *
1314                    (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1315            } else {
1316                int k, odind = cind + nci;
1317                if (uncst) {
1318                    int ione = 1, kk;
1319                    double *rr = Calloc(nci, double);
1320                    nlme_symmetrize(tmp, nci);
1321                    for (j = 0; j < nci; j++, cind++) {
1322                        for (k = 0; k < nci; k++) rr[k] = chol[j + k*nci];
1323                        REAL(ans)[cind] = 0.;
1324                        for (k = j; k < nci; k++) {
1325                            for (kk = j; kk < nci; kk++) {
1326                                REAL(ans)[cind] += rr[k] * rr[kk] *
1327                                    tmp[kk * nci + k];
1328                            }
1329                        }
1330                        for (k = 0; k < nci; k++) rr[k] *= rr[j];
1331                        for (k = j + 1; k < nci; k++) {
1332                            REAL(ans)[odind++] =
1333                                F77_CALL(ddot)(&nci, rr, &ione, tmp + k, &nci) +
1334                                F77_CALL(ddot)(&nci, rr, &ione,
1335                                               tmp + k*nci, &ione);
1336                        }
1337                    }
1338                    Free(rr);
1339                } else {
1340                    for (j = 0; j < nci; j++) {
1341                        REAL(ans)[cind++] = tmp[j * ncip1];
1342                        for (k = j + 1; k < nci; k++) {
1343                            REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1344                        }
1345                    }
1346                }
1347                cind = odind;
1348            }
1349            Free(tmp); Free(chol);
1350        }
1351      UNPROTECT(1);      UNPROTECT(1);
1352      return val;      return ans;
1353  }  }
1354    
   
1355  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)
1356  {  {
1357      SEXP val, b;      SEXP val, b;
# Line 1271  Line 1385 
1385              while (ff[j + nn] == lev) nn++;              while (ff[j + nn] == lev) nn++;
1386              F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,              F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
1387                              &one, mm + j, &nobs,                              &one, mm + j, &nobs,
1388                              bb + lev - 1, &nci,                              bb + (lev - 1) * nci, &nci,
1389                              &one, vv + j, &nobs);                              &one, vv + j, &nobs);
1390              j += nn;              j += nn;
1391          }          }
# Line 1279  Line 1393 
1393      UNPROTECT(2);      UNPROTECT(2);
1394      return val;      return val;
1395  }  }
1396    
1397    /**
1398     * Return the unscaled variances
1399     *
1400     * @param x pointer to an ssclme object
1401     *
1402     * @return a list similar to the Omega list with the unscaled variances
1403     */
1404    SEXP ssclme_variances(SEXP x)
1405    {
1406        SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1407        int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1408            i, nf = length(Omg);
1409    
1410        for (i = 0; i < nf; i++) {
1411            double *mm = REAL(VECTOR_ELT(Omg, i));
1412            int j, nci = nc[i];
1413    
1414            F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1415            if (j)                  /* shouldn't happen */
1416                error("DPOTRF returned error code %d on Omega[%d]",
1417                      j, i + 1);
1418            F77_CALL(dpotri)("U", &nci, mm, &nci, &j);
1419            if (j)                  /* shouldn't happen */
1420                error("DTRTRI returned error code %d on Omega[%d]",
1421                      j, i + 1);
1422            nlme_symmetrize(mm, nci);
1423        }
1424        UNPROTECT(1);
1425        return Omg;
1426    }
1427    
1428    SEXP ssclme_collapse(SEXP x)
1429    {
1430        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),
1431            Omega = GET_SLOT(x, Matrix_OmegaSym),
1432            Dim = GET_SLOT(x, Matrix_DimSym);
1433        int i, nf = length(Omega), nz = INTEGER(Dim)[1];
1434        SEXP copy[] = {Matrix_DSym, Matrix_DIsqrtSym, Matrix_DimSym,
1435                       Matrix_GpSym, Matrix_LIiSym, Matrix_LIpSym,
1436                       Matrix_LIxSym, Matrix_LiSym, Matrix_LpSym,
1437                       Matrix_LxSym, Matrix_OmegaSym, Matrix_ParentSym,
1438                       Matrix_LIxSym, Matrix_LiSym, Matrix_LpSym,
1439                       Matrix_bVarSym, Matrix_devianceSym,
1440                       Matrix_devCompSym, Matrix_iSym, Matrix_ncSym,
1441                       Matrix_statusSym, Matrix_pSym, Matrix_xSym};
1442    
1443        for (i = 0; i < 23; i++)
1444            SET_SLOT(ans, copy[i], duplicate(GET_SLOT(x, copy[i])));
1445    
1446        INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1447        SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));
1448        REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;
1449        SET_SLOT(ans, Matrix_RXXSym, allocMatrix(REALSXP, 1, 1));
1450        REAL(GET_SLOT(ans, Matrix_RXXSym))[0] = NA_REAL;
1451        SET_SLOT(ans, Matrix_ZtXSym, allocMatrix(REALSXP, nz, 1));
1452        SET_SLOT(ans, Matrix_RZXSym, allocMatrix(REALSXP, nz, 1));
1453        LOGICAL(GET_SLOT(ans, Matrix_statusSym))[0] = 0;
1454        UNPROTECT(1);
1455        return ans;
1456    }
1457    
1458    

Legend:
Removed from v.28  
changed lines
  Added in v.125

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