Diff of /pkg/src/ssclme.c

Parent Directory | Revision Log | Patch

	revision 244, Mon Jul 12 12:45:17 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 1528	Line 1635
1635	UNPROTECT(1);	UNPROTECT(1);
1636	return ans;	return ans;
1637	}	}
1638

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

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

---

root@r-forge.r-project.org	ViewVC Help
Powered by ViewVC 1.0.0