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