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