Add verbose option to omxComputeIterate
[openmx:openmx.git] / src / Compute.cpp
1 /*
2  *  Copyright 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 "omxDefines.h"
18 #include "Compute.h"
19 #include "omxState.h"
20 #include "omxExportBackendState.h"
21 #include "omxRFitFunction.h"
22
23 std::vector<int> FitContext::markMatrices;
24
25 void FitContext::init()
26 {
27         size_t numParam = varGroup->vars.size();
28         fit = 0;
29         est = new double[numParam];
30         grad = new double[numParam];
31         hess = new double[numParam * numParam];
32 }
33
34 FitContext::FitContext()
35 {
36         parent = NULL;
37         varGroup = Global->freeGroup[0];
38         init();
39
40         size_t numParam = varGroup->vars.size();
41         for (size_t v1=0; v1 < numParam; v1++) {
42                 est[v1] = Global->freeGroup[0]->vars[v1]->start;
43                 grad[v1] = nan("unset");
44                 for (size_t v2=0; v2 < numParam; v2++) {
45                         hess[v1 * numParam + v2] = nan("unset");
46                 }
47         }
48 }
49
50 FitContext::FitContext(FitContext *parent)
51 {
52         varGroup = parent->varGroup;
53         init();
54         fit = parent->fit;
55         size_t numParam = varGroup->vars.size();
56         memcpy(est, parent->est, sizeof(double) * numParam);
57         memcpy(grad, parent->grad, sizeof(double) * numParam);
58         memcpy(hess, parent->hess, sizeof(double) * numParam * numParam);
59 }
60
61 // arg to control what to copy? usually don't want everything TODO
62 FitContext::FitContext(FitContext *parent, FreeVarGroup *varGroup)
63 {
64         this->parent = parent;
65         this->varGroup = varGroup;
66         init();
67
68         FreeVarGroup *src = parent->varGroup;
69         FreeVarGroup *dest = varGroup;
70         size_t svars = parent->varGroup->vars.size();
71         size_t dvars = varGroup->vars.size();
72
73         size_t d1 = 0;
74         for (size_t s1=0; s1 < src->vars.size(); ++s1) {
75                 if (src->vars[s1] != dest->vars[d1]) continue;
76                 est[d1] = parent->est[s1];
77                 grad[d1] = parent->grad[s1];
78
79                 size_t d2 = 0;
80                 for (size_t s2=0; s2 < src->vars.size(); ++s2) {
81                         if (src->vars[s2] != dest->vars[d2]) continue;
82                         hess[d1 * dvars + d2] = parent->hess[s1 * svars + s2];
83                         if (++d2 == dvars) break;
84                 }
85
86                 if (++d1 == dvars) break;
87         }
88         if (d1 != dvars) error("Parent free parameter group is not a superset");
89
90         // pda(parent->est, 1, svars);
91         // pda(est, 1, dvars);
92         // pda(parent->grad, 1, svars);
93         // pda(grad, 1, dvars);
94         // pda(parent->hess, svars, svars);
95         // pda(hess, dvars, dvars);
96 }
97
98 void FitContext::copyParamToModel(omxMatrix *mat)
99 { copyParamToModel(mat->currentState); }
100
101 void FitContext::copyParamToModel(omxMatrix *mat, double *at)
102 { copyParamToModel(mat->currentState, at); }
103
104 void FitContext::updateParentAndFree()
105 {
106         FreeVarGroup *src = varGroup;
107         FreeVarGroup *dest = parent->varGroup;
108         size_t svars = varGroup->vars.size();
109         size_t dvars = parent->varGroup->vars.size();
110
111         size_t s1 = 0;
112         for (size_t d1=0; d1 < dest->vars.size(); ++d1) {
113                 if (dest->vars[d1] != src->vars[s1]) continue;
114                 parent->est[d1] = est[s1];
115                 parent->grad[d1] = grad[s1];
116
117                 size_t s2 = 0;
118                 for (size_t d2=0; d2 < dest->vars.size(); ++d2) {
119                         if (dest->vars[d2] != src->vars[s2]) continue;
120                         parent->hess[d1 * dvars + d2] = hess[s1 * svars + s2];
121                         if (++s2 == svars) break;
122                 }
123
124                 if (++s1 == svars) break;
125         }
126         
127         // pda(est, 1, svars);
128         // pda(parent->est, 1, dvars);
129         // pda(grad, 1, svars);
130         // pda(parent->grad, 1, dvars);
131         // pda(hess, svars, svars);
132         // pda(parent->hess, dvars, dvars);
133
134         delete this;
135 }
136
137 void FitContext::log(enum omxFFCompute what)
138 {
139         size_t count = varGroup->vars.size();
140         std::string buf;
141         if (what & FF_COMPUTE_FIT) buf += string_snprintf("fit: %.2g\n", fit);
142         if (what & FF_COMPUTE_ESTIMATE) {
143                 buf += "est: c(";
144                 for (size_t vx=0; vx < count; ++vx) {
145                         buf += string_snprintf("%.2g", est[vx]);
146                         if (vx < count - 1) buf += ", ";
147                 }
148                 buf += ")\n";
149         }
150         mxLogBig(buf);
151 }
152
153 void FitContext::cacheFreeVarDependencies()
154 {
155         omxState *os = globalState;
156         size_t numMats = os->matrixList.size();
157         size_t numAlgs = os->algebraList.size();
158
159         markMatrices.clear();
160         markMatrices.resize(numMats + numAlgs, 0);
161
162         // More efficient to use the appropriate group instead of the default group. TODO
163         FreeVarGroup *varGroup = Global->freeGroup[0];
164         for (size_t freeVarIndex = 0; freeVarIndex < varGroup->vars.size(); freeVarIndex++) {
165                 omxFreeVar* freeVar = varGroup->vars[freeVarIndex];
166                 int *deps   = freeVar->deps;
167                 int numDeps = freeVar->numDeps;
168                 for (int index = 0; index < numDeps; index++) {
169                         markMatrices[deps[index] + numMats] = 1;
170                 }
171         }
172 }
173
174 static void omxRepopulateRFitFunction(omxFitFunction* oo, double* x, int n)
175 {
176         omxRFitFunction* rFitFunction = (omxRFitFunction*)oo->argStruct;
177
178         SEXP theCall, estimate;
179
180         PROTECT(estimate = allocVector(REALSXP, n));
181         double *est = REAL(estimate);
182         for(int i = 0; i < n ; i++) {
183                 est[i] = x[i];
184         }
185
186         PROTECT(theCall = allocVector(LANGSXP, 4));
187
188         SETCAR(theCall, install("imxUpdateModelValues"));
189         SETCADR(theCall, rFitFunction->model);
190         SETCADDR(theCall, rFitFunction->flatModel);
191         SETCADDDR(theCall, estimate);
192
193         REPROTECT(rFitFunction->model = eval(theCall, R_GlobalEnv), rFitFunction->modelIndex);
194
195         UNPROTECT(2); // theCall, estimate
196 }
197
198 void FitContext::copyParamToModel(omxState* os)
199 {
200         copyParamToModel(os, est);
201 }
202
203 void FitContext::copyParamToModel(omxState* os, double *at)
204 {
205         size_t numParam = varGroup->vars.size();
206         if(OMX_DEBUG) {
207                 mxLog("Copying %d free parameter estimates to model %p", numParam, os);
208         }
209
210         if(numParam == 0) return;
211
212         size_t numMats = os->matrixList.size();
213         size_t numAlgs = os->algebraList.size();
214
215         os->computeCount++;
216
217         if(OMX_VERBOSE) {
218                 std::string buf;
219                 buf += string_snprintf("Call: %d.%d (%d) ", os->majorIteration, os->minorIteration, os->computeCount);
220                 buf += ("Estimates: [");
221                 for(size_t k = 0; k < numParam; k++) {
222                         buf += string_snprintf(" %f", at[k]);
223                 }
224                 buf += ("]\n");
225                 mxLogBig(buf);
226         }
227
228         for(size_t k = 0; k < numParam; k++) {
229                 omxFreeVar* freeVar = varGroup->vars[k];
230                 for(size_t l = 0; l < freeVar->locations.size(); l++) {
231                         omxFreeVarLocation *loc = &freeVar->locations[l];
232                         omxMatrix *matrix = os->matrixList[loc->matrix];
233                         int row = loc->row;
234                         int col = loc->col;
235                         omxSetMatrixElement(matrix, row, col, at[k]);
236                         if(OMX_DEBUG) {
237                                 mxLog("Setting location (%d, %d) of matrix %d to value %f for var %d",
238                                         row, col, loc->matrix, at[k], k);
239                         }
240                 }
241         }
242
243         if (RFitFunction) omxRepopulateRFitFunction(RFitFunction, at, numParam);
244
245         for(size_t i = 0; i < numMats; i++) {
246                 if (markMatrices[i]) {
247                         int offset = ~(i - numMats);
248                         omxMarkDirty(os->matrixList[offset]);
249                 }
250         }
251
252         for(size_t i = 0; i < numAlgs; i++) {
253                 if (markMatrices[i + numMats]) {
254                         omxMarkDirty(os->algebraList[i]);
255                 }
256         }
257
258         if (!os->childList) return;
259
260         for(int i = 0; i < Global->numChildren; i++) {
261                 copyParamToModel(os->childList[i]);
262         }
263 }
264
265 FitContext::~FitContext()
266 {
267         delete [] est;
268         delete [] grad;
269         delete [] hess;
270 }
271
272 omxFitFunction *FitContext::RFitFunction = NULL;
273
274 void FitContext::setRFitFunction(omxFitFunction *rff)
275 {
276         if (rff) {
277                 Global->numThreads = 1;
278                 if (RFitFunction) {
279                         error("You can only create 1 MxRFitFunction per independent model");
280                 }
281         }
282         RFitFunction = rff;
283 }
284
285 omxCompute::~omxCompute()
286 {}
287
288 void omxComputeOperation::initFromFrontend(SEXP rObj)
289 {
290         SEXP slotValue;
291         PROTECT(slotValue = GET_SLOT(rObj, install("free.group")));
292         int paramGroup = INTEGER(slotValue)[0];
293         varGroup = Global->freeGroup[paramGroup];
294 }
295
296 class omxComputeAssign : public omxComputeOperation {
297         typedef omxComputeOperation super;
298         std::vector< int > from;
299         std::vector< int > to;
300         std::vector< int > pmap;
301
302  public:
303         virtual void initFromFrontend(SEXP rObj);
304         virtual void compute(FitContext *fc);
305         virtual void reportResults(FitContext *fc, MxRList *out);
306 };
307
308 class omxComputeSequence : public omxCompute {
309         std::vector< omxCompute* > clist;
310
311  public:
312         virtual void initFromFrontend(SEXP rObj);
313         virtual void compute(FitContext *fc);
314         virtual void reportResults(FitContext *fc, MxRList *out);
315         virtual double getOptimizerStatus();
316         virtual ~omxComputeSequence();
317 };
318
319 class omxComputeIterate : public omxCompute {
320         std::vector< omxCompute* > clist;
321         int maxIter;
322         double tolerance;
323         bool verbose;
324
325  public:
326         virtual void initFromFrontend(SEXP rObj);
327         virtual void compute(FitContext *fc);
328         virtual void reportResults(FitContext *fc, MxRList *out);
329         virtual double getOptimizerStatus();
330         virtual ~omxComputeIterate();
331 };
332
333 class omxComputeOnce : public omxComputeOperation {
334         typedef omxComputeOperation super;
335         omxMatrix *fitMatrix;
336         omxExpectation *expectation;
337         const char *context;
338
339  public:
340         virtual void initFromFrontend(SEXP rObj);
341         virtual void compute(FitContext *fc);
342         virtual void reportResults(FitContext *fc, MxRList *out);
343 };
344
345 static class omxCompute *newComputeSequence()
346 { return new omxComputeSequence(); }
347
348 static class omxCompute *newComputeIterate()
349 { return new omxComputeIterate(); }
350
351 static class omxCompute *newComputeOnce()
352 { return new omxComputeOnce(); }
353
354 static class omxCompute *newComputeAssign()
355 { return new omxComputeAssign(); }
356
357 struct omxComputeTableEntry {
358         char name[32];
359         omxCompute *(*ctor)();
360 };
361
362 static const struct omxComputeTableEntry omxComputeTable[] = {
363         {"MxComputeEstimatedHessian", &newComputeEstimatedHessian},
364         {"MxComputeGradientDescent", &newComputeGradientDescent},
365         {"MxComputeSequence", &newComputeSequence },
366         {"MxComputeIterate", &newComputeIterate },
367         {"MxComputeOnce", &newComputeOnce },
368         {"MxComputeAssign", &newComputeAssign },
369 };
370
371 omxCompute *omxNewCompute(omxState* os, const char *type)
372 {
373         omxCompute *got = NULL;
374
375         for (size_t fx=0; fx < OMX_STATIC_ARRAY_SIZE(omxComputeTable); fx++) {
376                 const struct omxComputeTableEntry *entry = omxComputeTable + fx;
377                 if(strcmp(type, entry->name) == 0) {
378                         got = entry->ctor();
379                         break;
380                 }
381         }
382
383         if (!got) error("Compute %s is not implemented", type);
384
385         return got;
386 }
387
388 void omxComputeAssign::initFromFrontend(SEXP rObj)
389 {
390         super::initFromFrontend(rObj);
391
392         SEXP slotValue;
393
394         PROTECT(slotValue = GET_SLOT(rObj, install("from")));
395         for (int mx=0; mx < length(slotValue); ++mx) {
396                 from.push_back(INTEGER(slotValue)[mx]);
397         }
398
399         PROTECT(slotValue = GET_SLOT(rObj, install("to")));
400         for (int mx=0; mx < length(slotValue); ++mx) {
401                 to.push_back(INTEGER(slotValue)[mx]);
402         }
403
404         if (to.size() != from.size()) {
405                 error("omxComputeAssign: length of from and to must match");
406         }
407
408         pmap.assign(varGroup->vars.size(), -1);
409         size_t offset=0;
410         for (size_t mx=0; mx < from.size(); mx++) {
411                 int fmat = from[mx];
412                 int tmat = to[mx];
413                 for (size_t p1=0; p1 < varGroup->vars.size(); ++p1) {
414                         omxFreeVarLocation *l1 = varGroup->vars[p1]->getLocation(fmat);
415                         if (!l1) continue;
416                         bool found = FALSE;
417                         for (size_t r1=0; r1 < varGroup->vars.size(); ++r1) {
418                                 size_t p2 = (r1 + offset) % varGroup->vars.size();
419                                 omxFreeVarLocation *l2 = varGroup->vars[p2]->getLocation(tmat);
420                                 if (!l2) continue;
421                                 if (l1->row != l2->row || l1->col != l2->col) continue;
422                                 if (pmap[p1] != -1) {
423                                         error("omxComputeAssign: cannot copy %s to more than 1 place",
424                                               globalState->matrixList[~fmat]->name);
425                                 }
426                                 if (p1 == p2) error("omxComputeAssign: cannot copy %s to itself",
427                                                     globalState->matrixList[~fmat]->name);
428                                 pmap[p1] = p2;
429                                 offset = p2+1;  // probably a good guess
430                                 found = TRUE;
431                                 break;
432                         }
433                         if (!found) error("omxComputeAssign: %s[%d,%d] is fixed, cannot copy %s",
434                                           globalState->matrixList[~tmat]->name, 1+l1->row, 1+l1->col,
435                                           globalState->matrixList[~fmat]->name);
436                 }
437         }
438 }
439
440 void omxComputeAssign::compute(FitContext *fc)
441 {
442         for (size_t px=0; px < varGroup->vars.size(); ++px) {
443                 if (pmap[px] == -1) continue;
444                 fc->est[ pmap[px] ] = fc->est[px];
445         }
446         fc->copyParamToModel(globalState);
447 }
448
449 void omxComputeAssign::reportResults(FitContext *fc, MxRList *out)
450 {}
451
452 void omxComputeSequence::initFromFrontend(SEXP rObj)
453 {
454         SEXP slotValue;
455         PROTECT(slotValue = GET_SLOT(rObj, install("steps")));
456
457         for (int cx = 0; cx < length(slotValue); cx++) {
458                 SEXP step = VECTOR_ELT(slotValue, cx);
459                 SEXP s4class;
460                 PROTECT(s4class = STRING_ELT(getAttrib(step, install("class")), 0));
461                 omxCompute *compute = omxNewCompute(globalState, CHAR(s4class));
462                 compute->initFromFrontend(step);
463                 if (isErrorRaised(globalState)) break;
464                 clist.push_back(compute);
465         }
466 }
467
468 void omxComputeSequence::compute(FitContext *fc)
469 {
470         for (size_t cx=0; cx < clist.size(); ++cx) {
471                 FitContext *context = fc;
472                 if (fc->varGroup != clist[cx]->varGroup) {
473                         context = new FitContext(fc, clist[cx]->varGroup);
474                 }
475                 clist[cx]->compute(context);
476                 if (context != fc) context->updateParentAndFree();
477                 if (isErrorRaised(globalState)) break;
478         }
479 }
480
481 void omxComputeSequence::reportResults(FitContext *fc, MxRList *out)
482 {
483         for (size_t cx=0; cx < clist.size(); ++cx) {
484                 FitContext *context = fc;
485                 if (fc->varGroup != clist[cx]->varGroup) {
486                         context = new FitContext(fc, clist[cx]->varGroup);
487                 }
488                 clist[cx]->reportResults(context, out);
489                 if (context != fc) context->updateParentAndFree();
490                 if (isErrorRaised(globalState)) break;
491         }
492 }
493
494 double omxComputeSequence::getOptimizerStatus()
495 {
496         // for backward compatibility, not indended to work generally
497         for (size_t cx=0; cx < clist.size(); ++cx) {
498                 double got = clist[cx]->getOptimizerStatus();
499                 if (got != NA_REAL) return got;
500         }
501         return NA_REAL;
502 }
503
504 omxComputeSequence::~omxComputeSequence()
505 {
506         for (size_t cx=0; cx < clist.size(); ++cx) {
507                 delete clist[cx];
508         }
509 }
510
511 void omxComputeIterate::initFromFrontend(SEXP rObj)
512 {
513         SEXP slotValue;
514
515         PROTECT(slotValue = GET_SLOT(rObj, install("maxIter")));
516         maxIter = INTEGER(slotValue)[0];
517
518         PROTECT(slotValue = GET_SLOT(rObj, install("tolerance")));
519         tolerance = REAL(slotValue)[0];
520         if (tolerance <= 0) error("tolerance must be positive");
521
522         PROTECT(slotValue = GET_SLOT(rObj, install("steps")));
523
524         for (int cx = 0; cx < length(slotValue); cx++) {
525                 SEXP step = VECTOR_ELT(slotValue, cx);
526                 SEXP s4class;
527                 PROTECT(s4class = STRING_ELT(getAttrib(step, install("class")), 0));
528                 omxCompute *compute = omxNewCompute(globalState, CHAR(s4class));
529                 compute->initFromFrontend(step);
530                 if (isErrorRaised(globalState)) break;
531                 clist.push_back(compute);
532         }
533
534         PROTECT(slotValue = GET_SLOT(rObj, install("verbose")));
535         verbose = asLogical(slotValue);
536 }
537
538 void omxComputeIterate::compute(FitContext *fc)
539 {
540         int iter = 0;
541         double prevFit = 0;
542         double change = tolerance * 10;
543         while (1) {
544                 for (size_t cx=0; cx < clist.size(); ++cx) {
545                         FitContext *context = fc;
546                         if (fc->varGroup != clist[cx]->varGroup) {
547                                 context = new FitContext(fc, clist[cx]->varGroup);
548                         }
549                         clist[cx]->compute(context);
550                         if (context != fc) context->updateParentAndFree();
551                         if (isErrorRaised(globalState)) break;
552                 }
553                 if (prevFit != 0) {
554                         change = prevFit - fc->fit;
555                         if (verbose) mxLog("fit %.9g change %.9g", fc->fit, change);
556                 }
557                 prevFit = fc->fit;
558                 if (isErrorRaised(globalState) || ++iter > maxIter || fabs(change) < tolerance) break;
559         }
560 }
561
562 void omxComputeIterate::reportResults(FitContext *fc, MxRList *out)
563 {
564         for (size_t cx=0; cx < clist.size(); ++cx) {
565                 FitContext *context = fc;
566                 if (fc->varGroup != clist[cx]->varGroup) {
567                         context = new FitContext(fc, clist[cx]->varGroup);
568                 }
569                 clist[cx]->reportResults(context, out);
570                 if (context != fc) context->updateParentAndFree();
571                 if (isErrorRaised(globalState)) break;
572         }
573 }
574
575 double omxComputeIterate::getOptimizerStatus()
576 {
577         // for backward compatibility, not indended to work generally
578         for (size_t cx=0; cx < clist.size(); ++cx) {
579                 double got = clist[cx]->getOptimizerStatus();
580                 if (got != NA_REAL) return got;
581         }
582         return NA_REAL;
583 }
584
585 omxComputeIterate::~omxComputeIterate()
586 {
587         for (size_t cx=0; cx < clist.size(); ++cx) {
588                 delete clist[cx];
589         }
590 }
591
592 void omxComputeOnce::initFromFrontend(SEXP rObj)
593 {
594         super::initFromFrontend(rObj);
595
596         fitMatrix = omxNewMatrixFromSlot(rObj, globalState, "fitfunction");
597         if (fitMatrix) {
598                 setFreeVarGroup(fitMatrix->fitFunction, varGroup);
599                 omxCompleteFitFunction(fitMatrix);
600         }
601
602         SEXP slotValue;
603         PROTECT(slotValue = GET_SLOT(rObj, install("expectation")));
604         if (length(slotValue)) {
605                 int expNumber = INTEGER(slotValue)[0];  
606                 expectation = omxExpectationFromIndex(expNumber, globalState);
607                 setFreeVarGroup(expectation, varGroup);
608                 omxCompleteExpectation(expectation);
609         }
610
611         if (fitMatrix && expectation) {
612                 error("Cannot evaluate a fitfunction and expectation simultaneously");
613         }
614         if (!fitMatrix && !expectation) error("No function specified to evaluate");
615
616         PROTECT(slotValue = GET_SLOT(rObj, install("context")));
617         if (length(slotValue) == 0) {
618                 // OK
619         } else if (length(slotValue) == 1) {
620                 SEXP elem;
621                 PROTECT(elem = STRING_ELT(slotValue, 0));
622                 context = CHAR(elem);
623         }
624 }
625
626 void omxComputeOnce::compute(FitContext *fc)
627 {
628         if (fitMatrix) {
629                 for(size_t index = 0; index < globalState->matrixList.size(); index++) {
630                         omxMarkDirty(globalState->matrixList[index]);
631                 }
632                 for(size_t index = 0; index < globalState->algebraList.size(); index++) {
633                         omxMarkDirty(globalState->algebraList[index]);
634                 }
635                 omxFitFunctionCompute(fitMatrix->fitFunction, FF_COMPUTE_FIT, fc);
636                 fc->fit = fitMatrix->data[0];
637         } else if (expectation) {
638                 omxExpectationCompute(expectation, context);
639         }
640 }
641
642 void omxComputeOnce::reportResults(FitContext *fc, MxRList *out)
643 {
644         if (!fitMatrix) return;
645
646         omxPopulateFitFunction(fitMatrix, out);
647
648         out->push_back(std::make_pair(mkChar("minimum"), ScalarReal(fc->fit)));
649
650         size_t numFree = fc->varGroup->vars.size();
651         if (numFree) {
652                 SEXP estimate;
653                 PROTECT(estimate = allocVector(REALSXP, numFree));
654                 memcpy(REAL(estimate), fc->est, sizeof(double)*numFree);
655                 out->push_back(std::make_pair(mkChar("estimate"), estimate));
656         }
657 }