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