Removing unused variables that have been detected by gcc.
[openmx:openmx.git] / src / omxMLObjective.c
1  /*
2  *  Copyright 2007-2009 The OpenMx Project
3  *
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
7  *
8  *       http://www.apache.org/licenses/LICENSE-2.0
9  *
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.
15  */
16
17 #include <R.h>
18 #include <Rinternals.h>
19 #include <Rdefines.h>
20 #include <R_ext/Rdynload.h>
21 #include <R_ext/BLAS.h>
22 #include <R_ext/Lapack.h>
23 #include "omxAlgebraFunctions.h"
24
25 extern void omxInitFIMLObjective(omxObjective* oo, SEXP rObj); // Included in case ML gets a raw data object.
26 extern void omxCreateFIMLObjective(omxObjective* oo, SEXP rObj, omxMatrix* cov, omxMatrix* means);
27 void omxCreateMLObjective(omxObjective* oo, SEXP rObj, omxMatrix* cov, omxMatrix* means);
28
29 #ifndef _OMX_ML_OBJECTIVE_
30 #define _OMX_ML_OBJECTIVE_ TRUE
31
32 typedef struct omxMLObjective {
33
34         omxMatrix* observedCov;
35         omxMatrix* observedMeans;
36         omxMatrix* expectedCov;
37         omxMatrix* expectedMeans;
38         omxMatrix* localCov;
39         omxMatrix* localProd;
40         omxMatrix* P;
41         omxMatrix* C;
42         omxMatrix* I;
43         double n;
44         double logDetObserved;
45
46         double* work;
47         int lwork;
48
49 } omxMLObjective;
50
51 void omxDestroyMLObjective(omxObjective *oo) {
52
53         if(OMX_DEBUG) {Rprintf("Freeing ML Objective.");}
54         omxMLObjective* omlo = ((omxMLObjective*)oo->argStruct);
55
56         if(omlo->localCov != NULL)      omxFreeMatrixData(omlo->localCov);
57         if(omlo->localProd != NULL)     omxFreeMatrixData(omlo->localProd);
58         if(omlo->P != NULL)                     omxFreeMatrixData(omlo->P);
59         if(omlo->C != NULL)                     omxFreeMatrixData(omlo->C);
60         if(omlo->I != NULL)                     omxFreeMatrixData(omlo->I);
61 }
62
63 omxRListElement* omxSetFinalReturnsMLObjective(omxObjective *oo, int *numReturns) {
64         *numReturns = 3;
65         omxRListElement* retVal = (omxRListElement*) R_alloc(*numReturns, sizeof(omxRListElement));
66         double det = 0.0;
67         omxMatrix* cov = ((omxMLObjective*)oo->argStruct)->observedCov;
68         int ncols = ((omxMLObjective*)oo->argStruct)->observedCov->cols;
69     
70         retVal[0].numValues = 1;
71         retVal[0].values = (double*) R_alloc(1, sizeof(double));
72         strncpy(retVal[0].label, "Minus2LogLikelihood", 20);
73         retVal[0].values[0] = omxMatrixElement(oo->matrix, 0, 0);
74
75         retVal[1].numValues = 1;
76         retVal[1].values = (double*) R_alloc(1, sizeof(double));
77         strncpy(retVal[1].label, "SaturatedLikelihood", 20);
78         retVal[1].values[0] = (((omxMLObjective*)oo->argStruct)->logDetObserved + ncols) * (((omxMLObjective*)oo->argStruct)->n - 1);
79
80         retVal[2].numValues = 1;
81         retVal[2].values = (double*) R_alloc(1, sizeof(double));
82         strncpy(retVal[2].label, "IndependenceLikelihood", 23);
83         // Independence model assumes all-zero manifest covariances.
84         // (det(expected) + tr(observed * expected^-1)) * (n - 1);
85         // expected is the diagonal of the observed.  Inverse expected is 1/each diagonal value.
86         // Therefore the diagonal elements of observed * expected^-1 are each 1.
87         // So the trace of the matrix is the same as the number of columns.
88         // The determinant of a diagonal matrix is the product of the diagonal elements.
89         // Since these are the same in the expected as in the observed, we can get 'em direct.
90         for(int i = 0; i < ncols; i++) {
91                 // We sum logs instead of logging the product.
92                 det += log(omxMatrixElement(cov, i, i));
93         }
94         if(OMX_DEBUG) { Rprintf("det: %f, tr: %f, n= %d, total:%f\n", det, ncols, ((omxMLObjective*)oo->argStruct)->n, (ncols + det) * (((omxMLObjective*)oo->argStruct)->n - 1)); }
95         if(OMX_DEBUG) { omxPrint(cov, "Observed:"); }
96         retVal[2].values[0] = (ncols + det) * (((omxMLObjective*)oo->argStruct)->n - 1);
97
98         return retVal;
99 }
100
101 void omxCallMLObjective(omxObjective *oo) {     // TODO: Figure out how to give access to other per-iteration structures.
102
103         if(OMX_DEBUG) { Rprintf("Beginning ML Evaluation.\n");}
104         // Requires: Data, means, covariances.
105
106         double sum = 0.0, det = 0.0;
107         char u = 'U';
108         char r = 'R';
109         int info = 0;
110         double oned = 1.0;
111         double zerod = 0.0;
112         double minusoned = -1.0;
113         int onei = 1;
114         double fmean = 0.0;
115
116         omxMatrix *scov, *smeans, *cov, *means, *localCov, *localProd, *P, *C;
117
118         omxMLObjective *omo = ((omxMLObjective*)oo->argStruct);
119         
120     /* Locals for readability.  Compiler should cut through this. */
121         scov            = omo->observedCov;
122         smeans          = omo->observedMeans;
123         cov                     = omo->expectedCov;
124         means           = omo->expectedMeans;
125         localCov        = omo->localCov;
126         localProd       = omo->localProd;
127         P                       = omo->P;
128         C                       = omo->C;
129         double n        = omo->n;
130         double Q        = omo->logDetObserved;
131         omxObjective* subObjective = oo->subObjective;
132
133     /* Recompute and recopy */
134         if(!(subObjective == NULL)) {
135         omxObjectiveCompute(subObjective);
136         } else {
137                 omxRecompute(cov);
138                 omxRecompute(means);
139         }
140
141         omxCopyMatrix(localCov, cov);                           // But expected cov is destroyed in inversion
142
143         if(OMX_DEBUG_ALGEBRA) {
144                 omxPrint(scov, "Observed Covariance is");
145                 omxPrint(localCov, "Implied Covariance Is");
146                 omxPrint(cov, "Original Covariance Is");
147         }
148
149         /* Calculate |expected| */
150
151 //      F77_CALL(dgetrf)(&(localCov->cols), &(localCov->rows), localCov->data, &(localCov->cols), ipiv, &info);
152         F77_CALL(dpotrf)(&u, &(localCov->cols), localCov->data, &(localCov->cols), &info);
153
154         if(OMX_DEBUG_ALGEBRA) { Rprintf("Info on LU Decomp: %d\n", info);
155         omxPrint(localCov, "After Decomp:");}
156         if(info > 0) {
157                 char *errstr = calloc(250, sizeof(char));
158                 sprintf(errstr, "Expected covariance matrix is non-positive-definite");
159                 if(oo->matrix->currentState->computeCount <= 0) {
160                         strncat(errstr, " at starting values", 20);
161                 }
162                 strncat(errstr, ".\n", 3);
163                 omxRaiseError(oo->matrix->currentState, -1, errstr);                                            // Raise error
164                 free(errstr);
165                 return;                                                                                                                                         // Leave output untouched
166         }
167
168         //det = log(det)        // TVO: changed multiplication of det to first log and the summing up; this line should just set det to zero.
169         for(info = 0; info < localCov->cols; info++) {          // |cov| is the square of the product of the diagonal elements of U from the LU factorization.
170                 det += log(fabs(localCov->data[info+localCov->rows*info])); // TVO: changed * to + and added fabs command
171         }
172         det *= 2.0;             // TVO: instead of det *= det;
173
174         if(OMX_DEBUG_ALGEBRA) { Rprintf("Determinant of Expected Cov: %f\n", exp(det)); }
175         // TVO: removed det = log(fabs(det))
176         if(OMX_DEBUG_ALGEBRA) { Rprintf("Log of Determinant of Expected Cov: %f\n", det); }
177
178         /* Calculate Expected^(-1) */
179 //      F77_CALL(dgetri)(&(localCov->rows), localCov->data, &(localCov->cols), ipiv, work, lwork, &info);
180         F77_CALL(dpotri)(&u, &(localCov->rows), localCov->data, &(localCov->cols), &info);
181         if(OMX_DEBUG_ALGEBRA) { Rprintf("Info on Invert: %d\n", info); }
182
183         if(OMX_DEBUG_ALGEBRA) {omxPrint(cov, "Expected Covariance Matrix:");}
184         if(OMX_DEBUG_ALGEBRA) {omxPrint(localCov, "Inverted Matrix:");}
185
186         /* Calculate C = Observed * expected^(-1) */
187
188         if(OMX_DEBUG_ALGEBRA) {Rprintf("Call is: DSYMM(%d, %d, %f, %0x, %d, %0x, %d, %f, %0x, %d)",
189                                         (scov->rows), (localCov->cols), oned, scov->data, (localCov->leading),
190                                         localCov->data, (localCov->leading), zerod, localProd->data, (localProd->leading));}
191
192
193         // Stop gcc from issuing a warning
194         int majority = *(scov->majority) == 'n' ? scov->rows : scov->cols;
195
196         /*  TODO:  Make sure leading edges are being appropriately calculated, and sub them back into this */
197         F77_CALL(dsymm)(&r, &u, &(localCov->rows), &(scov->cols),
198                                         &oned, localCov->data, &(majority),
199                                         scov->data, &(majority),
200                                         &zerod, localProd->data, &(localProd->leading));
201
202     /* And get the trace of the result */
203
204         for(info = 0; info < localCov->cols; info++) {
205                 sum += localProd->data[info*localCov->cols + info];
206         }
207
208 //      for(info = 0; info < (localCov->cols * localCov->rows); info++) {
209 //              sum += localCov->data[info] * scov->data[info];
210 //      }
211
212         if(OMX_DEBUG_ALGEBRA) {omxPrint(scov, "Observed Covariance Matrix:");}
213         if(OMX_DEBUG_ALGEBRA) {omxPrint(localCov, "Inverse Matrix:");}
214         if(OMX_DEBUG_ALGEBRA) {omxPrint(localProd, "Product Matrix:");}
215
216         if(means != NULL) {
217                 if(OMX_DEBUG_ALGEBRA) { Rprintf("Means Likelihood Calculation"); }
218                 omxRecompute(means);
219                 omxCopyMatrix(P, means);
220                 // P = means - smeans
221                 if(OMX_DEBUG_ALGEBRA) {omxPrint(means, "means");}
222                 if(OMX_DEBUG_ALGEBRA) {omxPrint(smeans, "smeans");}
223                 F77_CALL(daxpy)(&(smeans->cols), &minusoned, smeans->data, &onei, P->data, &onei);
224                 if(OMX_DEBUG_ALGEBRA) {omxPrint(P, "means - smeans");}
225                 // C = P * Cov
226                 F77_CALL(dsymv)(&u, &(localCov->rows), &oned, localCov->data, &(localCov->leading), P->data, &onei, &zerod, C->data, &onei);
227                 // P = C * P'
228                 fmean = F77_CALL(ddot)(&(C->cols), P->data, &onei, C->data, &onei);
229
230                 if(OMX_DEBUG_ALGEBRA) { Rprintf("Mean contribution to likelihood is %f per row.\n", fmean); }
231                 if(fmean < 0.0) fmean = 0.0;
232         }
233
234         oo->matrix->data[0] = (sum + det) * (n - 1) + fmean * (n);
235
236         if(OMX_DEBUG) { Rprintf("MLObjective value comes to: %f (Chisq: %f).\n", oo->matrix->data[0], (sum + det) - Q - cov->cols); }
237
238 }
239
240 unsigned short int omxNeedsUpdateMLObjective(omxObjective* oo) {
241         if(omxNeedsUpdate(((omxMLObjective*)oo->argStruct)->expectedCov))
242                 return 1;
243         if(((omxMLObjective*)oo->argStruct)->expectedMeans != NULL)
244                 return omxNeedsUpdate(((omxMLObjective*)oo->argStruct)->expectedMeans);
245         return 0;
246 }
247
248 void omxPopulateMLAttributes(omxObjective *oo, SEXP algebra) {
249     if(OMX_DEBUG) { Rprintf("Populating ML Attributes.\n"); }
250
251         omxMLObjective *argStruct = ((omxMLObjective*)oo->argStruct);
252         omxMatrix *expCovInt = argStruct->expectedCov;                  // Expected covariance
253         omxMatrix *expMeanInt = argStruct->expectedMeans;                       // Expected means
254
255         SEXP expCovExt, expMeanExt;
256         PROTECT(expCovExt = allocMatrix(REALSXP, expCovInt->rows, expCovInt->cols));
257         for(int row = 0; row < expCovInt->rows; row++)
258                 for(int col = 0; col < expCovInt->cols; col++)
259                         REAL(expCovExt)[col * expCovInt->rows + row] =
260                                 omxMatrixElement(expCovInt, row, col);
261
262         if (expMeanInt != NULL) {
263                 PROTECT(expMeanExt = allocMatrix(REALSXP, expMeanInt->rows, expMeanInt->cols));
264                 for(int row = 0; row < expMeanInt->rows; row++)
265                         for(int col = 0; col < expMeanInt->cols; col++)
266                                 REAL(expMeanExt)[col * expMeanInt->rows + row] =
267                                         omxMatrixElement(expMeanInt, row, col);
268         } else {
269                 PROTECT(expMeanExt = allocMatrix(REALSXP, 0, 0));               
270         }   
271         setAttrib(algebra, install("expCov"), expCovExt);
272         setAttrib(algebra, install("expMean"), expMeanExt);
273         UNPROTECT(2);
274
275 }
276
277 void omxSetMLObjectiveCalls(omxObjective* oo) {
278         
279         /* Set Objective Calls to ML Objective Calls */
280         oo->objectiveFun = omxCallMLObjective;
281         oo->needsUpdateFun = omxNeedsUpdateMLObjective;
282         oo->destructFun = omxDestroyMLObjective;
283         oo->setFinalReturns = omxSetFinalReturnsMLObjective;
284         oo->populateAttrFun = omxPopulateMLAttributes;
285         
286 }
287
288 void omxInitMLObjective(omxObjective* oo, SEXP rObj) {
289
290         if(OMX_DEBUG) { Rprintf("Initializing ML objective function.\n"); }
291
292         SEXP nextMatrix;
293         omxMatrix *cov, *means;
294         
295         /* Read and set expected means and variances */
296         PROTECT(nextMatrix = GET_SLOT(rObj, install("means")));
297         if(OMX_DEBUG) { Rprintf("Processing Expected Means.\n"); }
298         if(!R_FINITE(INTEGER(nextMatrix)[0])) {
299                 if(OMX_DEBUG) {
300                         Rprintf("ML: No Expected Means.\n");
301                 }
302                 means = NULL;
303         } else {
304                 means = omxNewMatrixFromMxIndex(nextMatrix, oo->matrix->currentState);
305                 if(OMX_DEBUG) { Rprintf("Means matrix created at 0x%x.\n", means); }
306         }
307         UNPROTECT(1);
308
309         PROTECT(nextMatrix = GET_SLOT(rObj, install("covariance")));
310         if(OMX_DEBUG) { Rprintf("Processing Expected Covariance.\n"); }
311         cov = omxNewMatrixFromMxIndex(nextMatrix, oo->matrix->currentState);
312         UNPROTECT(1);
313         
314         omxCreateMLObjective(oo, rObj, cov, means);
315         
316 }
317
318 void omxCreateMLObjective(omxObjective* oo, SEXP rObj, omxMatrix* cov, omxMatrix* means) {
319         
320         SEXP nextMatrix;
321         int info = 0;
322         double det = 1.0;
323         char u = 'U';
324         
325         omxSetMLObjectiveCalls(oo);
326         
327         if(OMX_DEBUG) { Rprintf("Retrieving data.\n"); }
328         PROTECT(nextMatrix = GET_SLOT(rObj, install("data")));
329         omxData* dataMat = omxNewDataFromMxDataPtr(nextMatrix, oo->matrix->currentState);
330         if(strncmp(omxDataType(dataMat), "cov", 3) != 0 && strncmp(omxDataType(dataMat), "cor", 3) != 0) {
331                 if(strncmp(omxDataType(dataMat), "raw", 3) == 0) {
332                         if(OMX_DEBUG) { Rprintf("Raw Data: Converting to FIML.\n"); }
333             UNPROTECT(1); // Cleanup before sending to omxCreateFIML.
334                         omxCreateFIMLObjective(oo, rObj, cov, means);
335                         return;
336                 }
337                 char *errstr = calloc(250, sizeof(char));
338                 sprintf(errstr, "ML Objective unable to handle data type %s.\n", omxDataType(dataMat));
339                 omxRaiseError(oo->matrix->currentState, -1, errstr);
340                 free(errstr);
341                 if(OMX_DEBUG) { Rprintf("ML Objective unable to handle data type %s.  Aborting.\n", omxDataType(dataMat)); }
342                 return;
343         }
344         
345         if(cov == NULL) {
346                 omxRaiseError(oo->matrix->currentState, OMX_DEVELOPER_ERROR,
347                         "Developer Error in ML-based objective object: ML-based subobjectives must specify a model-implied covariance matrix.\nIf you are not developing a new objective type, you should probably post this to the OpenMx forums.");
348                 return;
349         }
350
351         omxMLObjective *newObj = (omxMLObjective*) R_alloc(1, sizeof(omxMLObjective));
352
353         newObj->expectedCov = cov;
354         newObj->expectedMeans = means;
355
356         if(OMX_DEBUG) { Rprintf("Processing Observed Covariance.\n"); }
357         newObj->observedCov = omxDataMatrix(dataMat, NULL);
358         if(OMX_DEBUG) { Rprintf("Processing Observed Means.\n"); }
359         newObj->observedMeans = omxDataMeans(dataMat, NULL, NULL);
360         if(OMX_DEBUG && newObj->observedMeans == NULL) { Rprintf("ML: No Observed Means.\n"); }
361         if(OMX_DEBUG) { Rprintf("Processing n.\n"); }
362         newObj->n = omxDataNumObs(dataMat);
363         UNPROTECT(1); // nextMatrix
364         
365         // Error Checking: Observed/Expected means must agree.  
366         // ^ is XOR: true when one is false and the other is not.
367         if((newObj->expectedMeans == NULL) ^ (newObj->observedMeans == NULL)) {
368             if(newObj->expectedMeans != NULL) {
369                     omxRaiseError(oo->matrix->currentState, OMX_ERROR,
370                             "Observed means not detected, but an expected means matrix was specified.\n  If you provide observed means, you must specify a model for the means.\n");
371                     return;
372             } else {
373                     omxRaiseError(oo->matrix->currentState, OMX_ERROR,
374                             "Observed means were provided, but an expected means matrix was not specified.\n  If you  wish to model the means, you must provide observed means.\n");
375                     return;             
376             }
377         }
378
379         /* Temporary storage for calculation */
380         int rows = newObj->observedCov->rows;
381         int cols = newObj->observedCov->cols;
382         newObj->localCov = omxInitMatrix(NULL, rows, cols, TRUE, oo->matrix->currentState);
383         newObj->localProd = omxInitMatrix(NULL, rows, cols, TRUE, oo->matrix->currentState);
384         newObj->P = omxInitMatrix(NULL, 1, cols, TRUE, oo->matrix->currentState);
385         newObj->C = omxInitMatrix(NULL, rows, cols, TRUE, oo->matrix->currentState);
386         newObj->I = omxInitMatrix(NULL, rows, cols, TRUE, oo->matrix->currentState);
387
388         for(int i = 0; i < rows; i++) omxSetMatrixElement(newObj->I, i, i, 1.0);
389
390         omxCopyMatrix(newObj->localCov, newObj->observedCov);
391
392         newObj->lwork = newObj->expectedCov->rows;
393         newObj->work = (double*)R_alloc(newObj->lwork, sizeof(double));
394
395         F77_CALL(dpotrf)(&u, &(newObj->localCov->cols), newObj->localCov->data, &(newObj->localCov->cols), &info);
396
397         if(OMX_DEBUG) { Rprintf("Info on LU Decomp: %d\n", info); }
398         if(info != 0) {
399                 char *errstr = calloc(250, sizeof(char));
400                 sprintf(errstr, "Observed Covariance Matrix is non-positive-definite.\n");
401                 omxRaiseError(oo->matrix->currentState, -1, errstr);
402                 free(errstr);
403                 return;
404         }
405         for(info = 0; info < newObj->localCov->cols; info++) {
406                 det *= omxMatrixElement(newObj->localCov, info, info);
407         }
408         det *= det;                                     // Product of squares.
409
410         if(OMX_DEBUG) { Rprintf("Determinant of Observed Cov: %f\n", det); }
411         newObj->logDetObserved = log(det);
412         if(OMX_DEBUG) { Rprintf("Log Determinant of Observed Cov: %f\n", newObj->logDetObserved); }
413
414         omxCopyMatrix(newObj->localCov, newObj->expectedCov);
415     oo->argStruct = (void*)newObj;
416
417 }
418
419 #endif /* _OMX_ML_OBJECTIVE_ */