Sanity check shape of FitFunction@expectation slot
[openmx:openmx.git] / src / omxFitFunction.c
1 /*
2  *  Copyright 2007-2012 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 /***********************************************************
18
19 *  omxFitFunction.cc
20 *
21 *  Created: Timothy R. Brick    Date: 2008-11-13 12:33:06
22 *
23 *       FitFunction objects are a subclass of data matrix that evaluates
24 *   itself anew at each iteration, so that any changes to
25 *   free parameters can be incorporated into the update.
26 *   // Question: Should FitFunction be a ``subtype'' of 
27 *   // omxAlgebra or a separate beast entirely?
28 *
29 **********************************************************/
30
31 #include "omxFitFunction.h"
32
33 typedef struct omxFitFunctionTableEntry omxFitFunctionTableEntry;
34
35 struct omxFitFunctionTableEntry {
36
37         char name[32];
38         void (*initFun)(omxFitFunction*, SEXP);
39
40 };
41
42 extern void omxInitAlgebraFitFunction(omxFitFunction *off, SEXP rObj);
43 extern void omxInitWLSFitFunction(omxFitFunction *off, SEXP rObj);
44 extern void omxInitRowFitFunction(omxFitFunction *off, SEXP rObj);
45 extern void omxInitMLFitFunction(omxFitFunction *off, SEXP rObj);
46 extern void omxInitRFitFunction(omxFitFunction *off, SEXP rObj);
47
48 static const omxFitFunctionTableEntry omxFitFunctionSymbolTable[] = {
49         {"MxFitFunctionAlgebra",                        &omxInitAlgebraFitFunction},
50         {"MxFitFunctionWLS",                            &omxInitWLSFitFunction},
51         {"MxFitFunctionRow",                            &omxInitRowFitFunction},
52         {"MxFitFunctionML",                             &omxInitMLFitFunction},
53         {"MxFitFunctionR",                                      &omxInitRFitFunction},
54         {"", 0}
55 };
56
57 void omxCalculateStdErrorFromHessian(double scale, omxFitFunction *off) {
58         /* This function calculates the standard errors from the hessian matrix */
59         // sqrt(diag(solve(hessian)))
60
61         if(off->hessian == NULL) return;
62         
63         int numParams = off->matrix->currentState->numFreeParams;
64         
65         if(off->stdError == NULL) {
66                 off->stdError = (double*) R_alloc(numParams, sizeof(double));
67         }
68         
69         double* stdErr = off->stdError;
70         
71         double* hessian = off->hessian;
72         double* workspace = (double *) Calloc(numParams * numParams, double);
73         
74         for(int i = 0; i < numParams; i++)
75                 for(int j = 0; j <= i; j++)
76                         workspace[i*numParams+j] = hessian[i*numParams+j];              // Populate upper triangle
77         
78         char u = 'U';
79         int ipiv[numParams];
80         int lwork = -1;
81         double temp;
82         int info = 0;
83         
84         F77_CALL(dsytrf)(&u, &numParams, workspace, &numParams, ipiv, &temp, &lwork, &info);
85         
86         lwork = (temp > numParams?temp:numParams);
87         
88         double* work = (double*) Calloc(lwork, double);
89         
90         F77_CALL(dsytrf)(&u, &numParams, workspace, &numParams, ipiv, work, &lwork, &info);
91         
92         if(info != 0) {
93                 
94                 off->stdError = NULL;
95                 
96         } else {
97                 
98                 F77_CALL(dsytri)(&u, &numParams, workspace, &numParams, ipiv, work, &info);
99         
100                 if(info != 0) {
101                         off->stdError = NULL;
102                 } else {
103                         for(int i = 0; i < numParams; i++) {
104                                 stdErr[i] = sqrt(scale) * sqrt(workspace[i * numParams + i]);
105                         }
106                 }
107         }
108         
109         Free(workspace);
110         Free(work);
111         
112 }
113
114 void omxInitEmptyFitFunction(omxFitFunction *off) {
115         /* Sets everything to NULL to avoid bad pointer calls */
116         
117         memset(off, 0, sizeof(omxFitFunction));
118 }
119
120 void omxFreeFitFunctionArgs(omxFitFunction *off) {
121         if(off==NULL) return;
122     
123         /* Completely destroy the fit function structures */
124         if(OMX_DEBUG) {Rprintf("Freeing fit function object at 0x%x.\n", off);}
125         if(off->matrix != NULL) {
126                 if(off->destructFun != NULL) {
127                         if(OMX_DEBUG) {Rprintf("Calling fit function destructor for 0x%x.\n", off);}
128                         off->destructFun(off);
129                 }
130                 off->matrix = NULL;
131         }
132 }
133
134 void omxFitFunctionCompute(omxFitFunction *off) {
135         if(OMX_DEBUG_ALGEBRA) { 
136             Rprintf("FitFunction compute: 0x%0x (needed: %s).\n", off, (off->matrix->isDirty?"Yes":"No"));
137         }
138
139         off->computeFun(off);
140
141         omxMarkClean(off->matrix);
142
143 }
144
145 void omxDuplicateFitMatrix(omxMatrix *tgt, const omxMatrix *src, omxState* newState) {
146
147         if(tgt == NULL || src == NULL) return;
148         if(src->fitFunction == NULL) return;
149     
150         omxFillMatrixFromMxFitFunction(tgt, src->fitFunction->rObj, src->hasMatrixNumber, src->matrixNumber);
151
152 }
153
154 omxFitFunction* omxCreateDuplicateFitFunction(omxFitFunction *tgt, const omxFitFunction *src, omxState* newState) {
155
156         if(OMX_DEBUG) {Rprintf("Duplicating fit function 0x%x into 0x%x.", src, tgt);}
157
158         if(src == NULL) {
159                 return NULL;
160         }
161         
162         if(tgt == NULL) {
163         tgt = (omxFitFunction*) R_alloc(1, sizeof(omxFitFunction));
164         omxInitEmptyFitFunction(tgt);
165     } else {
166                 omxRaiseError(newState, -1,
167                         "omxCreateDuplicateFitFunction requested to overwrite target");
168                 return NULL;
169         }
170
171         memcpy(tgt, src, sizeof(omxFitFunction));
172         return tgt;
173
174 }
175
176 void omxFillMatrixFromMxFitFunction(omxMatrix* om, SEXP rObj,
177         unsigned short hasMatrixNumber, int matrixNumber) {
178
179         SEXP slotValue, fitFunctionClass;
180         omxFitFunction *obj = (omxFitFunction*) R_alloc(1, sizeof(omxFitFunction));
181         omxInitEmptyFitFunction(obj);
182
183         /* Register FitFunction and Matrix with each other */
184         obj->matrix = om;
185         omxResizeMatrix(om, 1, 1, FALSE);                                       // FitFunction matrices MUST be 1x1.
186         om->fitFunction = obj;
187         om->hasMatrixNumber = hasMatrixNumber;
188         om->matrixNumber = matrixNumber;
189         
190         /* Get FitFunction Type */
191         PROTECT(fitFunctionClass = STRING_ELT(getAttrib(rObj, install("class")), 0));
192         {
193           const char *fitType = CHAR(fitFunctionClass);
194         
195           /* Switch based on fit function type. */ 
196           const omxFitFunctionTableEntry *entry = omxFitFunctionSymbolTable;
197           while (entry->initFun) {
198             if(strncmp(fitType, entry->name, MAX_STRING_LEN) == 0) {
199               obj->fitType = entry->name;
200               obj->initFun = entry->initFun;
201               break;
202             }
203             entry += 1;
204           }
205
206           if(obj->initFun == NULL) {
207             const int MaxErrorLen = 256;
208             char newError[MaxErrorLen];
209             snprintf(newError, MaxErrorLen, "Fit function %s not implemented.\n", fitType);
210             omxRaiseError(om->currentState, -1, newError);
211             return;
212           }
213         }
214         UNPROTECT(1);   /* fitType */
215
216         PROTECT(slotValue = GET_SLOT(rObj, install("expectation")));
217         if (LENGTH(slotValue) != 1) {
218             const int MaxErrorLen = 256;
219             char newError[MaxErrorLen];
220             snprintf(newError, MaxErrorLen, "Fit function %s expectation improperly initialized\n", obj->fitType);
221             error(newError);
222         }
223         int expNumber = INTEGER(slotValue)[0];  
224         if(expNumber == NA_INTEGER) {                                           // Has no expectation associated with it
225                 obj->expectation = NULL;
226         } else {
227                 obj->expectation = omxNewExpectationFromExpectationIndex(expNumber, om->currentState);
228         }
229         UNPROTECT(1);   /* slotValue */
230         
231         obj->rObj = rObj;
232         obj->initFun(obj, rObj);
233
234         if(obj->computeFun == NULL) {// If initialization fails, error code goes in argStruct
235                 char *errorCode;
236                 if(om->currentState->statusCode != 0) {
237                         errorCode = om->currentState->statusMsg;
238                 } else {
239                         // If no error code is reported, we report that.
240                         errorCode = "No error code reported.";
241                 }
242                 if(obj->argStruct != NULL) {
243                         errorCode = (char*)(obj->argStruct);
244                 }
245         const int MaxErrorLen = 256;
246         char newError[MaxErrorLen];
247         snprintf(newError, MaxErrorLen, "Could not initialize fit function %s.  Error: %s\n",
248                         obj->fitType, errorCode); 
249                 omxRaiseError(om->currentState, -1, newError);
250         }
251         
252         obj->matrix->isDirty = TRUE;
253
254 }
255
256 void omxFitFunctionGradient(omxFitFunction* off, double* gradient) {
257         if(!(off->gradientFun == NULL)) { off->gradientFun(off, gradient); }
258         return;
259 }
260
261 void omxFitFunctionPrint(omxFitFunction* off, char* d) {
262         Rprintf("(FitFunction, type %s) ", off->fitType);
263         omxPrintMatrix(off->matrix, d);
264 }
265
266
267 /* Helper functions */
268 omxMatrix* omxNewMatrixFromIndexSlot(SEXP rObj, omxState* currentState, char* const slotName) {
269         SEXP slotValue;
270         omxMatrix* newMatrix = NULL;
271         if(strncmp(slotName, "", 1) == 0) return NULL;
272         PROTECT(slotValue = GET_SLOT(rObj, install(slotName)));
273         newMatrix = omxNewMatrixFromMxIndex(slotValue, currentState);
274         if(newMatrix != NULL) omxRecompute(newMatrix);
275         else if(OMX_DEBUG) Rprintf("No slot %s found.\n", slotName);
276         UNPROTECT(1);
277         return newMatrix;
278 }
279
280 omxData* omxNewDataFromDataSlot(SEXP rObj, omxState* currentState, char* const dataSlotName) {
281         
282         SEXP slotValue;
283         
284         PROTECT(slotValue = GET_SLOT(rObj, install(dataSlotName)));
285         if(OMX_DEBUG) { Rprintf("Data Element %d.\n", AS_INTEGER(slotValue)); }
286         omxData* dataElt = omxNewDataFromMxDataPtr(slotValue, currentState);
287         UNPROTECT(1); // newMatrix
288         
289         return dataElt;
290         
291 }