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