2 * Copyright 2007-2013 The OpenMx Project
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
8 * http://www.apache.org/licenses/LICENSE-2.0
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
18 #include <Rinternals.h>
20 #include <R_ext/Rdynload.h>
21 #include <R_ext/BLAS.h>
22 #include <R_ext/Lapack.h>
23 #include "omxAlgebraFunctions.h"
24 #include "omxWLSFitFunction.h"
26 void flattenDataToVector(omxMatrix* cov, omxMatrix* means, omxThresholdColumn* thresholds, int nThresholds, omxMatrix* vector) {
27 // TODO: vectorize data flattening
29 for(int j = 0; j < cov->rows; j++) {
30 for(int k = 0; k <= j; k++) {
31 omxSetVectorElement(vector, nextLoc, omxMatrixElement(cov, k, j)); // Use upper triangle in case of SYMM-style mat.
36 for(int j = 0; j < cov->rows; j++) {
37 omxSetVectorElement(vector, nextLoc, omxVectorElement(means, j));
41 if (thresholds != NULL) {
42 for(int j = 0; j < nThresholds; j++) {
43 omxThresholdColumn* thresh = thresholds + j;
44 for(int k = 0; k < thresh->numThresholds; k++) {
45 omxSetVectorElement(vector, nextLoc, omxMatrixElement(thresh->matrix, k, thresh->column));
52 void omxDestroyWLSFitFunction(omxFitFunction *oo) {
54 if(OMX_DEBUG) {mxLog("Freeing WLS FitFunction.");}
55 omxWLSFitFunction* owo = ((omxWLSFitFunction*)oo->argStruct);
57 if(owo->observedFlattened != NULL) omxFreeMatrixData(owo->observedFlattened);
58 if(owo->expectedFlattened != NULL) omxFreeMatrixData(owo->expectedFlattened);
59 if(owo->weights != NULL) omxFreeMatrixData(owo->weights);
60 if(owo->B != NULL) omxFreeMatrixData(owo->B);
61 if(owo->P != NULL) omxFreeMatrixData(owo->P);
64 omxRListElement* omxSetFinalReturnsWLSFitFunction(omxFitFunction *oo, int *numReturns) {
66 omxRListElement* retVal = (omxRListElement*) R_alloc(*numReturns, sizeof(omxRListElement));
68 retVal[0].numValues = 1;
69 retVal[0].values = (double*) R_alloc(1, sizeof(double));
70 strncpy(retVal[0].label, "Minus2LogLikelihood", 20);
71 retVal[0].values[0] = omxMatrixElement(oo->matrix, 0, 0);
73 retVal[1].numValues = 1;
74 retVal[1].values = (double*) R_alloc(1, sizeof(double));
75 strncpy(retVal[1].label, "SaturatedLikelihood", 20);
76 retVal[1].values[0] = 0;
78 retVal[2].numValues = 1;
79 retVal[2].values = (double*) R_alloc(1, sizeof(double));
80 strncpy(retVal[2].label, "IndependenceLikelihood", 23);
81 retVal[2].values[0] = 0;
83 retVal[3].numValues = 1;
84 retVal[3].values = (double*) R_alloc(1, sizeof(double));
85 strncpy(retVal[3].label, "ADFMisfit", 20);
86 retVal[3].values[0] = omxMatrixElement(oo->matrix, 0, 0);
88 // retVal[1].numValues = 1;
89 // retVal[1].values = (double*) R_alloc(1, sizeof(double));
90 // strncpy(retVal[1].label, "SaturatedADFMisfit", 20);
91 // retVal[1].values[0] = (((omxWLSFitFunction*)oo->argStruct)->logDetObserved + ncols) * (((omxWLSFitFunction*)oo->argStruct)->n - 1);
93 // retVal[2].numValues = 1;
94 // retVal[2].values = (double*) R_alloc(1, sizeof(double));
95 // strncpy(retVal[2].label, "IndependenceLikelihood", 23);
96 // // Independence model assumes all-zero manifest covariances.
97 // // (det(expected) + tr(observed * expected^-1)) * (n - 1);
98 // // expected is the diagonal of the observed. Inverse expected is 1/each diagonal value.
99 // // Therefore the diagonal elements of observed * expected^-1 are each 1.
100 // // So the trace of the matrix is the same as the number of columns.
101 // // The determinant of a diagonal matrix is the product of the diagonal elements.
102 // // Since these are the same in the expected as in the observed, we can get 'em direct.
103 // for(int i = 0; i < ncols; i++) {
104 // // We sum logs instead of logging the product.
105 // det += log(omxMatrixElement(cov, i, i));
107 // if(OMX_DEBUG) { mxLog("det: %f, tr: %f, n= %d, total:%f", det, ncols, ((omxWLSFitFunction*)oo->argStruct)->n, (ncols + det) * (((omxWLSFitFunction*)oo->argStruct)->n - 1)); }
108 // if(OMX_DEBUG) { omxPrint(cov, "Observed:"); }
109 // retVal[2].values[0] = (ncols + det) * (((omxWLSFitFunction*)oo->argStruct)->n - 1);
114 static void omxCallWLSFitFunction(omxFitFunction *oo, int want, double *gradient, double *hessian) { // TODO: Figure out how to give access to other per-iteration structures.
116 if(OMX_DEBUG) { mxLog("Beginning WLS Evaluation.");}
117 // Requires: Data, means, covariances.
121 omxMatrix *oCov, *oMeans, *eCov, *eMeans, *P, *B, *weights, *oFlat, *eFlat;
123 omxThresholdColumn *oThresh, *eThresh;
125 omxWLSFitFunction *owo = ((omxWLSFitFunction*)oo->argStruct);
127 /* Locals for readability. Compiler should cut through this. */
128 oCov = owo->observedCov;
129 oMeans = owo->observedMeans;
130 oThresh = owo->observedThresholds;
131 eCov = owo->expectedCov;
132 eMeans = owo->expectedMeans;
133 eThresh = owo->expectedThresholds;
134 oFlat = owo->observedFlattened;
135 eFlat = owo->expectedFlattened;
136 weights = owo->weights;
139 int nThresh = owo->nThresholds;
142 omxExpectation* expectation = oo->expectation;
144 /* Recompute and recopy */
145 omxExpectationCompute(expectation);
147 flattenDataToVector(oCov, oMeans, oThresh, nThresh, oFlat);
148 flattenDataToVector(eCov, eMeans, eThresh, nThresh, eFlat);
150 omxCopyMatrix(B, oFlat);
152 omxDAXPY(-1.0, eFlat, B);
154 if(weights != NULL) {
155 omxDGEMV(TRUE, 1.0, weights, B, 0.0, P);
157 // ULS Case: Memcpy faster than dgemv.
158 // TODO: Better to use an omxMatrix duplicator here.
159 memcpy(P, B, B->cols*sizeof(double));
161 sum = F77_CALL(ddot)(&(P->cols), P->data, &onei, B->data, &onei);
163 oo->matrix->data[0] = sum;
165 if(OMX_DEBUG) { mxLog("WLSFitFunction value comes to: %f.", oo->matrix->data[0]); }
169 void omxPopulateWLSAttributes(omxFitFunction *oo, SEXP algebra) {
170 if(OMX_DEBUG) { mxLog("Populating WLS Attributes."); }
172 omxWLSFitFunction *argStruct = ((omxWLSFitFunction*)oo->argStruct);
173 omxMatrix *expCovInt = argStruct->expectedCov; // Expected covariance
174 omxMatrix *expMeanInt = argStruct->expectedMeans; // Expected means
175 omxMatrix *weightInt = argStruct->weights; // Expected means
177 SEXP expCovExt, expMeanExt, weightExt, gradients;
178 PROTECT(expCovExt = allocMatrix(REALSXP, expCovInt->rows, expCovInt->cols));
179 for(int row = 0; row < expCovInt->rows; row++)
180 for(int col = 0; col < expCovInt->cols; col++)
181 REAL(expCovExt)[col * expCovInt->rows + row] =
182 omxMatrixElement(expCovInt, row, col);
184 if (expMeanInt != NULL) {
185 PROTECT(expMeanExt = allocMatrix(REALSXP, expMeanInt->rows, expMeanInt->cols));
186 for(int row = 0; row < expMeanInt->rows; row++)
187 for(int col = 0; col < expMeanInt->cols; col++)
188 REAL(expMeanExt)[col * expMeanInt->rows + row] =
189 omxMatrixElement(expMeanInt, row, col);
191 PROTECT(expMeanExt = allocMatrix(REALSXP, 0, 0));
194 PROTECT(weightExt = allocMatrix(REALSXP, weightInt->rows, weightInt->cols));
195 for(int row = 0; row < weightInt->rows; row++)
196 for(int col = 0; col < weightInt->cols; col++)
197 REAL(weightExt)[col * weightInt->rows + row] =
198 omxMatrixElement(weightInt, row, col);
201 if(0) { /* TODO fix for new internal API
202 int nLocs = Global->numFreeParams;
203 double gradient[Global->numFreeParams];
204 for(int loc = 0; loc < nLocs; loc++) {
205 gradient[loc] = NA_REAL;
207 //oo->gradientFun(oo, gradient);
208 PROTECT(gradients = allocMatrix(REALSXP, 1, nLocs));
210 for(int loc = 0; loc < nLocs; loc++)
211 REAL(gradients)[loc] = gradient[loc];
214 PROTECT(gradients = allocMatrix(REALSXP, 0, 0));
217 setAttrib(algebra, install("expCov"), expCovExt);
218 setAttrib(algebra, install("expMean"), expMeanExt);
219 setAttrib(algebra, install("weights"), weightExt);
220 setAttrib(algebra, install("gradients"), gradients);
226 void omxSetWLSFitFunctionCalls(omxFitFunction* oo) {
228 /* Set FitFunction Calls to WLS FitFunction Calls */
229 oo->computeFun = omxCallWLSFitFunction;
230 oo->destructFun = omxDestroyWLSFitFunction;
231 oo->setFinalReturns = omxSetFinalReturnsWLSFitFunction;
232 oo->populateAttrFun = omxPopulateWLSAttributes;
235 void omxInitWLSFitFunction(omxFitFunction* oo) {
237 omxMatrix *cov, *means, *weights;
239 if(OMX_DEBUG) { mxLog("Initializing WLS FitFunction function."); }
243 omxSetWLSFitFunctionCalls(oo);
245 if(OMX_DEBUG) { mxLog("Retrieving expectation.\n"); }
246 if (!oo->expectation) { error("%s requires an expectation", oo->fitType); }
248 if(OMX_DEBUG) { mxLog("Retrieving data.\n"); }
249 omxData* dataMat = oo->expectation->data;
251 if(strncmp(omxDataType(dataMat), "acov", 4) != 0 && strncmp(omxDataType(dataMat), "cov", 3) != 0) {
252 char *errstr = (char*) calloc(250, sizeof(char));
253 sprintf(errstr, "WLS FitFunction unable to handle data type %s. Data must be of type 'acov'.\n", omxDataType(dataMat));
254 omxRaiseError(oo->matrix->currentState, -1, errstr);
256 if(OMX_DEBUG) { mxLog("WLS FitFunction unable to handle data type %s. Aborting.", omxDataType(dataMat)); }
260 omxWLSFitFunction *newObj = (omxWLSFitFunction*) R_alloc(1, sizeof(omxWLSFitFunction));
262 /* Get Expectation Elements */
263 newObj->expectedCov = omxGetExpectationComponent(oo->expectation, oo, "cov");
264 newObj->expectedMeans = omxGetExpectationComponent(oo->expectation, oo, "means");
265 newObj->nThresholds = oo->expectation->numOrdinal;
266 newObj->expectedThresholds = oo->expectation->thresholds;
267 // FIXME: threshold structure should be asked for by omxGetExpectationComponent
269 /* Read and set expected means, variances, and weights */
270 cov = omxDataMatrix(dataMat, NULL);
271 means = omxDataMeans(dataMat, NULL, NULL);
272 weights = omxDataAcov(dataMat, NULL);
274 newObj->observedCov = cov;
275 newObj->observedMeans = means;
276 newObj->weights = weights;
277 newObj->n = omxDataNumObs(dataMat);
278 newObj->nThresholds = omxDataNumFactor(dataMat);
281 // Error Checking: Observed/Expected means must agree.
282 // ^ is XOR: true when one is false and the other is not.
283 if((newObj->expectedMeans == NULL) ^ (newObj->observedMeans == NULL)) {
284 if(newObj->expectedMeans != NULL) {
285 omxRaiseError(oo->matrix->currentState, OMX_ERROR,
286 "Observed means not detected, but an expected means matrix was specified.\n If you wish to model the means, you must provide observed means.\n");
289 omxRaiseError(oo->matrix->currentState, OMX_ERROR,
290 "Observed means were provided, but an expected means matrix was not specified.\n If you provide observed means, you must specify a model for the means.\n");
295 if((newObj->expectedThresholds == NULL) ^ (newObj->observedThresholds == NULL)) {
296 if(newObj->expectedMeans != NULL) {
297 omxRaiseError(oo->matrix->currentState, OMX_ERROR,
298 "Observed thresholds not detected, but an expected thresholds matrix was specified.\n If you wish to model the thresholds, you must provide observed thresholds.\n ");
301 omxRaiseError(oo->matrix->currentState, OMX_ERROR,
302 "Observed thresholds were provided, but an expected thresholds matrix was not specified.\nIf you provide observed thresholds, you must specify a model for the thresholds.\n");
307 /* Error check weight matrix size */
308 int ncol = newObj->observedCov->cols;
309 vectorSize = (ncol * (ncol + 1) ) / 2;
310 if(newObj->expectedMeans != NULL) {
311 vectorSize = vectorSize + ncol;
313 if(newObj->observedThresholds != NULL) {
314 for(int i = 0; i < newObj->nThresholds; i++) {
315 vectorSize = vectorSize + newObj->observedThresholds[i].numThresholds;
319 if(weights != NULL && weights->rows != weights->cols && weights->cols != vectorSize) {
320 omxRaiseError(oo->matrix->currentState, OMX_DEVELOPER_ERROR,
321 "Developer Error in WLS-based FitFunction object: WLS-based expectation specified an incorrectly-sized weight matrix.\nIf you are not developing a new expectation type, you should probably post this to the OpenMx forums.");
326 // FIXME: More error checking for incoming Fit Functions
328 /* Temporary storage for calculation */
329 newObj->observedFlattened = omxInitMatrix(NULL, vectorSize, 1, TRUE, oo->matrix->currentState);
330 newObj->expectedFlattened = omxInitMatrix(NULL, vectorSize, 1, TRUE, oo->matrix->currentState);
331 newObj->P = omxInitMatrix(NULL, 1, vectorSize, TRUE, oo->matrix->currentState);
332 newObj->B = omxInitMatrix(NULL, vectorSize, 1, TRUE, oo->matrix->currentState);
334 flattenDataToVector(newObj->observedCov, newObj->observedMeans, newObj->observedThresholds, newObj->nThresholds, newObj->observedFlattened);
335 flattenDataToVector(newObj->expectedCov, newObj->expectedMeans, newObj->expectedThresholds, newObj->nThresholds, newObj->expectedFlattened);
337 oo->argStruct = (void*)newObj;