Store algebraList in std::vector
[openmx:openmx.git] / src / omxState.cpp
1 /*
2  *  Copyright 2007-2013 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 *  omxState.cc
20 *
21 *  Created: Timothy R. Brick    Date: 2009-06-05
22 *
23 *       omxStates carry the current optimization state
24 *
25 **********************************************************/
26
27 #include <stdarg.h>
28 #include "omxState.h"
29
30 /* Initialize and Destroy */
31         void omxInitState(omxState* state, omxState *parentState) {
32                 state->ciMaxIterations = 5;
33                 state->numThreads = 1;
34                 state->numHessians = 0;
35                 state->calculateStdErrors = FALSE;
36
37                 state->numExpects = 0;
38                 state->numConstraints = 0;
39                 state->numFreeParams = 0;
40                 state->numChildren = 0;
41                 state->childList = NULL;
42                 state->expectationList = NULL;
43                 state->parentState = parentState;
44                 state->conList = NULL;
45                 state->freeVarList = NULL;
46
47                 state->majorIteration = 0;
48                 state->minorIteration = 0;
49                 state->startTime = 0;
50                 state->endTime = 0;
51                 state->numCheckpoints = 0;
52                 state->checkpointList = NULL;
53                 state->chkptText1 = NULL;
54                 state->chkptText2 = NULL;
55
56                 state->computeCount = -1;
57                 state->currentRow = -1;
58
59                 strncpy(state->statusMsg, "", 1);
60         }
61
62         omxState* omxGetState(omxState* os, int stateNumber) {
63                 // TODO: Need to implement a smarter way to enumerate children
64                 if(stateNumber == 0) return os;
65                 if((stateNumber-1) < os->numChildren) {
66                         return(os->childList[stateNumber-1]);
67                 } else {
68                         error("Not implemented");
69                         // TODO: Account for unequal numbers of grandchild states
70                         int subState = (stateNumber - os->numChildren - 1);
71                         return omxGetState(os->childList[subState % os->numChildren], subState / os->numChildren);
72                 }
73         }
74
75         void omxSetMajorIteration(omxState *state, int value) {
76                 state->majorIteration = value;
77                 for(int i = 0; i < state->numChildren; i++) {
78                         omxSetMajorIteration(state->childList[i], value);
79                 }
80         }
81
82         void omxSetMinorIteration(omxState *state, int value) {
83                 state->minorIteration = value;
84                 for(int i = 0; i < state->numChildren; i++) {
85                         omxSetMinorIteration(state->childList[i], value);
86                 }
87         }
88         
89         void omxDuplicateState(omxState* tgt, omxState* src) {
90                 tgt->numExpects                 = src->numExpects;
91                 tgt->dataList                   = src->dataList;
92                 tgt->numChildren                = 0;
93                 
94                 // Duplicate matrices and algebras and build parentLists.
95                 tgt->parentState                = src;
96                 tgt->expectationList    = (omxExpectation**) R_alloc(tgt->numExpects, sizeof(omxExpectation*));
97                 tgt->markMatrices               = src->markMatrices; // TODO, unused in children?
98                                 
99                 memset(tgt->expectationList, 0, sizeof(omxExpectation*) * tgt->numExpects);
100
101                 for(size_t mx = 0; mx < src->matrixList.size(); mx++) {
102                         // TODO: Smarter inference for which matrices to duplicate
103                         tgt->matrixList.push_back(omxDuplicateMatrix(src->matrixList[mx], tgt));
104                 }
105
106                 tgt->numConstraints     = src->numConstraints;
107                 tgt->conList                    = (omxConstraint*) R_alloc(tgt->numConstraints, sizeof(omxConstraint));
108                 for(int j = 0; j < tgt->numConstraints; j++) {
109                         tgt->conList[j].size   = src->conList[j].size;
110                         tgt->conList[j].opCode = src->conList[j].opCode;
111                         tgt->conList[j].lbound = src->conList[j].lbound;
112                         tgt->conList[j].ubound = src->conList[j].ubound;
113                         tgt->conList[j].result = omxDuplicateMatrix(src->conList[j].result, tgt);
114                 }
115
116                 for(size_t j = 0; j < src->algebraList.size(); j++) {
117                         // TODO: Smarter inference for which algebras to duplicate
118                         tgt->algebraList.push_back(omxDuplicateMatrix(src->algebraList[j], tgt));
119                 }
120
121                 for(int j = 0; j < tgt->numExpects; j++) {
122                         // TODO: Smarter inference for which expectations to duplicate
123                         tgt->expectationList[j] = omxDuplicateExpectation(src->expectationList[j], tgt);
124                 }
125
126                 for(size_t j = 0; j < tgt->algebraList.size(); j++) {
127                         omxDuplicateAlgebra(tgt->algebraList[j], src->algebraList[j], tgt);
128                 }
129
130                 for(int j = 0; j < tgt->numExpects; j++) {
131                         // TODO: Smarter inference for which expectations to duplicate
132                         omxCompleteExpectation(tgt->expectationList[j]);
133                 }
134
135                 tgt->childList                  = NULL;
136
137                 tgt->numFreeParams                      = src->numFreeParams;
138                 tgt->freeVarList                = new omxFreeVar[tgt->numFreeParams];
139                 for(int j = 0; j < tgt->numFreeParams; j++) {
140                         int numDeps                                                     = src->freeVarList[j].numDeps;
141
142                         tgt->freeVarList[j].lbound                      = src->freeVarList[j].lbound;
143                         tgt->freeVarList[j].ubound                      = src->freeVarList[j].ubound;
144                         tgt->freeVarList[j].locations                   = src->freeVarList[j].locations;
145                         tgt->freeVarList[j].numDeps                     = numDeps;
146                         
147                         tgt->freeVarList[j].deps                        = (int*) R_alloc(numDeps, sizeof(int));
148
149                         tgt->freeVarList[j].name                = src->freeVarList[j].name;
150
151                         for(int k = 0; k < numDeps; k++) {
152                                 tgt->freeVarList[j].deps[k] = src->freeVarList[j].deps[k];
153                         }
154                 }
155                 
156                 tgt->majorIteration     = 0;
157                 tgt->minorIteration     = 0;
158                 tgt->startTime                  = src->startTime;
159                 tgt->endTime                    = 0;
160                 
161                 // TODO: adjust checkpointing based on parallelization method
162                 tgt->numCheckpoints     = 0;
163                 tgt->checkpointList     = NULL;
164                 tgt->chkptText1                 = NULL;
165                 tgt->chkptText2                 = NULL;
166                                   
167                 tgt->computeCount               = src->computeCount;
168                 tgt->currentRow                 = src->currentRow;
169
170                 strncpy(tgt->statusMsg, "", 1);
171         }
172
173         omxMatrix* omxLookupDuplicateElement(omxState* os, omxMatrix* element) {
174                 if(os == NULL || element == NULL) return NULL;
175
176                 if (element->hasMatrixNumber) {
177                         int matrixNumber = element->matrixNumber;
178                         if (matrixNumber >= 0) {
179                                 return(os->algebraList[matrixNumber]);
180                         } else {
181                                 return(os->matrixList[-matrixNumber - 1]);
182                         }
183                 }
184
185                 omxConstraint* parentConList = os->parentState->conList;
186
187                 for(int i = 0; i < os->numConstraints; i++) {
188                         if(parentConList[i].result == element) {
189                                 if(os->conList[i].result != NULL) {   // Not sure of proper failure behavior here.
190                 return(os->conList[i].result);
191                                 } else {
192                     omxRaiseError(os, -2, "Initialization Copy Error: Constraint required but not yet processed.");
193             }
194                         }
195                 }
196
197                 return NULL;
198         }
199         
200         omxExpectation* omxLookupDuplicateExpectation(omxState* os, omxExpectation* ox) {
201                 if(os == NULL || ox == NULL) return NULL;
202
203                 return(os->expectationList[ox->expNum]);
204         }
205
206         void omxFreeState(omxState *state) {
207                 int k;
208
209                 if (state->numChildren > 0) {
210                         for(k = 0; k < state->numChildren; k++) {
211                                 // Data are not modified and not copied. The same memory
212                                 // is shared across all instances of state. We only need
213                                 // to free the data once, so let the parent do it.
214                                 state->childList[k]->dataList.clear();
215
216                                 omxFreeState(state->childList[k]);
217                         }
218                         Free(state->childList);
219                         state->childList = NULL;
220                         state->numChildren = 0;
221                 }
222
223                 for(size_t ax = 0; ax < state->algebraList.size(); ax++) {
224                         if(OMX_DEBUG) { Rprintf("Freeing Algebra %d at 0x%x.\n", ax, state->algebraList[ax]); }
225                         omxFreeAllMatrixData(state->algebraList[ax]);
226                 }
227
228                 if(OMX_DEBUG) { Rprintf("Freeing %d Matrices.\n", state->matrixList.size());}
229                 for(size_t mk = 0; mk < state->matrixList.size(); mk++) {
230                         if(OMX_DEBUG) { Rprintf("Freeing Matrix %d at 0x%x.\n", mk, state->matrixList[mk]); }
231                         omxFreeAllMatrixData(state->matrixList[mk]);
232                 }
233                 
234                 if(OMX_DEBUG) { Rprintf("Freeing %d Model Expectations.\n", state->numExpects);}
235                 for(k = 0; k < state->numExpects; k++) {
236                         if(OMX_DEBUG) { Rprintf("Freeing Expectation %d at 0x%x.\n", k, state->expectationList[k]); }
237                         omxFreeExpectationArgs(state->expectationList[k]);
238                 }
239
240                 if(OMX_DEBUG) { Rprintf("Freeing %d Constraints.\n", state->numConstraints);}
241                 for(k = 0; k < state->numConstraints; k++) {
242                         if(OMX_DEBUG) { Rprintf("Freeing Constraint %d at 0x%x.\n", k, state->conList[k]); }
243                         omxFreeAllMatrixData(state->conList[k].result);
244                 }
245
246                 if(OMX_DEBUG) { Rprintf("Freeing %d Data Sets.\n", state->dataList.size());}
247                 for(size_t dx = 0; dx < state->dataList.size(); dx++) {
248                         if(OMX_DEBUG) { Rprintf("Freeing Data Set %d at 0x%x.\n", dx, state->dataList[dx]); }
249                         omxFreeData(state->dataList[dx]);
250                 }
251
252                 delete [] state->freeVarList;
253
254         if(OMX_DEBUG) {Rprintf("Freeing %d Children.\n", state->numChildren);}
255         for(k = 0; k < state->numChildren; k++) {
256                         if(OMX_DEBUG) { Rprintf("Freeing Child State %d at 0x%x.\n", k, state->childList[k]); }
257                         omxFreeState(state->childList[k]);            
258         }
259
260                 if(OMX_DEBUG) { Rprintf("Freeing %d Checkpoints.\n", state->numCheckpoints);}
261                 for(k = 0; k < state->numCheckpoints; k++) {
262                         if(OMX_DEBUG) { Rprintf("Freeing Data Set %d at 0x%x.\n", k, state->checkpointList[k]); }
263                         omxCheckpoint oC = state->checkpointList[k];
264                         switch(oC.type) {
265                                 case OMX_FILE_CHECKPOINT:
266                                         fclose(oC.file);
267                                         break;
268                                 case OMX_CONNECTION_CHECKPOINT: // NYI :::DEBUG:::
269                                         // Do nothing: this should be handled by R upon return.
270                                         break;
271                         }
272                         if(state->chkptText1 != NULL) {
273                                 Free(state->chkptText1);
274                         }
275                         if(state->chkptText2 != NULL) {
276                                 Free(state->chkptText2);
277                         }
278                         // Checkpoint list itself is freed by R.
279                 }
280
281                 delete state;
282
283                 if(OMX_DEBUG) { Rprintf("State Freed.\n");}
284         }
285
286         void omxResetStatus(omxState *state) {
287                 int numChildren = state->numChildren;
288                 state->statusMsg[0] = '\0';
289                 for(int i = 0; i < numChildren; i++) {
290                         omxResetStatus(state->childList[i]);
291                 }
292         }
293
294 void omxRaiseErrorf(omxState *state, const char* errorMsg, ...)
295 {
296         va_list ap;
297         va_start(ap, errorMsg);
298         int fit = vsnprintf(state->statusMsg, MAX_STRING_LEN, errorMsg, ap);
299         va_end(ap);
300         if(OMX_DEBUG) {
301                 if (!(fit > -1 && fit < MAX_STRING_LEN)) {
302                         Rprintf("Error exceeded maximum length: %s\n", errorMsg);
303                 } else {
304                         Rprintf("Error raised: %s\n", state->statusMsg);
305                 }
306         }
307 }
308
309         void omxRaiseError(omxState *state, int errorCode, const char* errorMsg) { // DEPRECATED
310                 if(OMX_DEBUG && errorCode) { Rprintf("Error %d raised: %s\n", errorCode, errorMsg);}
311                 if(OMX_DEBUG && !errorCode) { Rprintf("Error status cleared."); }
312                 strncpy(state->statusMsg, errorMsg, 249);
313                 state->statusMsg[249] = '\0';
314         }
315
316         void omxStateNextRow(omxState *state) {
317                 state->currentRow++;
318         };
319
320         void omxStateNextEvaluation(omxState *state) {
321                 state->currentRow = -1;
322                 state->computeCount++;
323         };
324
325         void omxWriteCheckpointHeader(omxState *os, omxCheckpoint* oC) {
326                 // FIXME: Is it faster to allocate this on the stack?
327                 os->chkptText1 = (char*) Calloc((24 + 15 * os->numFreeParams), char);
328                 os->chkptText2 = (char*) Calloc(1.0 + 15.0 * os->numFreeParams*
329                         (os->numFreeParams + 1.0) / 2.0, char);
330                 if (oC->type == OMX_FILE_CHECKPOINT) {
331                         fprintf(oC->file, "iterations\ttimestamp\tobjective\t");
332                         for(int j = 0; j < os->numFreeParams; j++) {
333                                 if(strcmp(os->freeVarList[j].name, CHAR(NA_STRING)) == 0) {
334                                         fprintf(oC->file, "%s", os->freeVarList[j].name);
335                                 } else {
336                                         fprintf(oC->file, "\"%s\"", os->freeVarList[j].name);
337                                 }
338                                 if (j != os->numFreeParams - 1) fprintf(oC->file, "\t");
339                         }
340                         fprintf(oC->file, "\n");
341                         fflush(oC->file);
342                 }
343         }
344  
345         void omxWriteCheckpointMessage(omxState *os, char *msg) {
346                 for(int i = 0; i < os->numCheckpoints; i++) {
347                         omxCheckpoint* oC = &(os->checkpointList[i]);
348                         if(os->chkptText1 == NULL) {    // First one: set up output
349                                 omxWriteCheckpointHeader(os, oC);
350                         }
351                         if (oC->type == OMX_FILE_CHECKPOINT) {
352                                 fprintf(oC->file, "%d \"%s\" NA ", os->majorIteration, msg);
353                                 for(int j = 0; j < os->numFreeParams; j++) {
354                                         fprintf(oC->file, "NA ");
355                                 }
356                                 fprintf(oC->file, "\n");
357                         }
358                 }
359         }
360
361         void omxSaveCheckpoint(omxState *os, double* x, double* f, int force) {
362                 time_t now = time(NULL);
363                 int soFar = now - os->startTime;                // Translated into minutes
364                 int n;
365                 for(int i = 0; i < os->numCheckpoints; i++) {
366                         n = 0;
367                         omxCheckpoint* oC = &(os->checkpointList[i]);
368                         // Check based on time            
369                         if((oC->time > 0 && (soFar - oC->lastCheckpoint) >= oC->time) || force) {
370                                 oC->lastCheckpoint = soFar;
371                                 n = 1;
372                         }
373                         // Or iterations
374                         if((oC->numIterations > 0 && (os->majorIteration - oC->lastCheckpoint) >= oC->numIterations) || force) {
375                                 oC->lastCheckpoint = os->majorIteration;
376                                 n = 1;
377                         }
378
379                         if(n) {         //In either case, save a checkpoint.
380                                 if(os->chkptText1 == NULL) {    // First one: set up output
381                                         omxWriteCheckpointHeader(os, oC);
382                                 }
383                                 char tempstring[25];
384                                 sprintf(tempstring, "%d", os->majorIteration);
385
386                                 if(strncmp(os->chkptText1, tempstring, strlen(tempstring))) {   // Returns zero if they're the same.
387                                         struct tm * nowTime = localtime(&now);                                          // So this only happens if the text is out of date.
388                                         strftime(tempstring, 25, "%b %d %Y %I:%M:%S %p", nowTime);
389                                         sprintf(os->chkptText1, "%d \"%s\" %9.5f", os->majorIteration, tempstring, f[0]);
390                                         for(int j = 0; j < os->numFreeParams; j++) {
391                                                 sprintf(tempstring, " %9.5f", x[j]);
392                                                 strncat(os->chkptText1, tempstring, 14);
393                                         }
394                                 }
395
396                                 if(oC->type == OMX_FILE_CHECKPOINT) {
397                                         fprintf(oC->file, "%s", os->chkptText1);
398                                         if(oC->saveHessian)
399                                                 fprintf(oC->file, "%s", os->chkptText2);
400                                         fprintf(oC->file, "\n");
401                                         fflush(oC->file);
402                                 } else if(oC->type == OMX_CONNECTION_CHECKPOINT) {
403                                         warning("NYI: R_connections are not yet implemented.");
404                                         oC->numIterations = 0;
405                                         oC->time = 0;
406                                 }
407                         }
408                 }
409         }
410
411 void omxExamineFitOutput(omxState *state, omxMatrix *fitMatrix, int *mode)
412 {
413         if (!R_FINITE(fitMatrix->data[0])) {
414                 omxRaiseErrorf(state, "Fit function returned %g at iteration %d.%d",
415                                fitMatrix->data[0], state->majorIteration, state->minorIteration);
416                 *mode = -1;
417         }
418 }