Generic ComputeEM with Ramsay (1975) acceleration
[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 #include "matrix.h"
23 #include "omxBuffer.h"
24
25 void FitContext::init()
26 {
27         size_t numParam = varGroup->vars.size();
28         wanted = 0;
29         mac = parent? parent->mac : 0;
30         fit = parent? parent->fit : 0;
31         caution = parent? parent->caution : 0;
32         est = new double[numParam];
33         flavor = new int[numParam];
34         grad = new double[numParam];
35         hess = new double[numParam * numParam];
36         ihess = new double[numParam * numParam];
37         changedEstimates = false;
38         inform = INFORM_UNINITIALIZED;
39         iterations = 0;
40 }
41
42 FitContext::FitContext(std::vector<double> &startingValues)
43 {
44         parent = NULL;
45         varGroup = Global->freeGroup[0];
46         init();
47
48         size_t numParam = varGroup->vars.size();
49         if (startingValues.size() != numParam) {
50                 error("Got %d starting values for %d parameters",
51                       startingValues.size(), numParam);
52         }
53         memcpy(est, startingValues.data(), sizeof(double) * numParam);
54
55         for (size_t v1=0; v1 < numParam; v1++) {
56                 grad[v1] = nan("unset");
57                 for (size_t v2=0; v2 < numParam; v2++) {
58                         hess[v1 * numParam + v2] = nan("unset");
59                 }
60         }
61 }
62
63 FitContext::FitContext(FitContext *parent, FreeVarGroup *varGroup)
64 {
65         this->parent = parent;
66         this->varGroup = varGroup;
67         init();
68
69         FreeVarGroup *src = parent->varGroup;
70         FreeVarGroup *dest = varGroup;
71         size_t svars = parent->varGroup->vars.size();
72         size_t dvars = varGroup->vars.size();
73         if (dvars == 0) return;
74         mapToParent.resize(dvars);
75
76         size_t d1 = 0;
77         for (size_t s1=0; s1 < src->vars.size(); ++s1) {
78                 if (src->vars[s1] != dest->vars[d1]) continue;
79                 mapToParent[d1] = s1;
80                 est[d1] = parent->est[s1];
81
82                 if (parent->wanted & (FF_COMPUTE_GRADIENT | FF_COMPUTE_HESSIAN | FF_COMPUTE_INFO)) {
83                         grad[d1] = parent->grad[s1];
84
85                         size_t d2 = 0;
86                         for (size_t s2=0; s2 < src->vars.size(); ++s2) {
87                                 if (src->vars[s2] != dest->vars[d2]) continue;
88                                 hess[d1 * dvars + d2] = parent->hess[s1 * svars + s2];
89                                 if (++d2 == dvars) break;
90                         }
91                 }
92
93                 // ihess TODO?
94
95                 if (++d1 == dvars) break;
96         }
97         if (d1 != dvars) error("Parent free parameter group is not a superset");
98
99         // pda(parent->est, 1, svars);
100         // pda(est, 1, dvars);
101         // pda(parent->grad, 1, svars);
102         // pda(grad, 1, dvars);
103         // pda(parent->hess, svars, svars);
104         // pda(hess, dvars, dvars);
105 }
106
107 void FitContext::copyParamToModel(omxMatrix *mat)
108 { copyParamToModel(mat->currentState); }
109
110 void FitContext::copyParamToModel(omxMatrix *mat, double *at)
111 { copyParamToModel(mat->currentState, at); }
112
113 void FitContext::updateParent()
114 {
115         FreeVarGroup *src = varGroup;
116         FreeVarGroup *dest = parent->varGroup;
117         size_t svars = varGroup->vars.size();
118         size_t dvars = parent->varGroup->vars.size();
119
120         parent->wanted |= wanted;
121         parent->fit = fit;
122         parent->mac = mac;
123         parent->caution = caution;
124
125         // rewrite using mapToParent TODO
126
127         if (svars > 0) {
128                 size_t s1 = 0;
129                 for (size_t d1=0; d1 < dest->vars.size(); ++d1) {
130                         if (dest->vars[d1] != src->vars[s1]) continue;
131                         parent->est[d1] = est[s1];
132
133                         if (wanted & (FF_COMPUTE_GRADIENT | FF_COMPUTE_HESSIAN | FF_COMPUTE_INFO)) {
134                                 parent->grad[d1] = grad[s1];
135
136                                 size_t s2 = 0;
137                                 for (size_t d2=0; d2 < dest->vars.size(); ++d2) {
138                                         if (dest->vars[d2] != src->vars[s2]) continue;
139                                         parent->hess[d1 * dvars + d2] = hess[s1 * svars + s2];
140                                         if (++s2 == svars) break;
141                                 }
142                         }
143
144                         // ihess TODO?
145
146                         if (++s1 == svars) break;
147                 }
148                 if (wanted & FF_COMPUTE_PARAMFLAVOR) {
149                         for (size_t s1=0; s1 < src->vars.size(); ++s1) {
150                                 parent->flavor[mapToParent[s1]] = flavor[s1];
151                         }
152                 }
153         }
154         
155         // pda(est, 1, svars);
156         // pda(parent->est, 1, dvars);
157         // pda(grad, 1, svars);
158         // pda(parent->grad, 1, dvars);
159         // pda(hess, svars, svars);
160         // pda(parent->hess, dvars, dvars);
161 }
162
163 void FitContext::updateParentAndFree()
164 {
165         updateParent();
166         delete this;
167 }
168
169 void FitContext::log(const char *where)
170 {
171         log(where, wanted);
172 }
173
174 void FitContext::log(const char *where, int what)
175 {
176         size_t count = varGroup->vars.size();
177         std::string buf(where);
178         buf += " ---\n";
179         if (what & FF_COMPUTE_MAXABSCHANGE) buf += string_snprintf("MAC: %.5f\n", mac);
180         if (what & FF_COMPUTE_FIT) buf += string_snprintf("fit: %.5f\n", fit);
181         if (what & FF_COMPUTE_ESTIMATE) {
182                 buf += string_snprintf("est %lu: c(", count);
183                 for (size_t vx=0; vx < count; ++vx) {
184                         buf += string_snprintf("%.5f", est[vx]);
185                         if (vx < count - 1) buf += ", ";
186                 }
187                 buf += ")\n";
188         }
189         if (what & FF_COMPUTE_GRADIENT) {
190                 buf += string_snprintf("grad %lu: c(", count);
191                 for (size_t vx=0; vx < count; ++vx) {
192                         buf += string_snprintf("%.5f", grad[vx]);
193                         if (vx < count - 1) buf += ", ";
194                 }
195                 buf += ")\n";
196         }
197         if (what & (FF_COMPUTE_HESSIAN | FF_COMPUTE_INFO)) {
198                 buf += string_snprintf("hess %lux%lu: c(", count, count);
199                 for (size_t v1=0; v1 < count; ++v1) {
200                         for (size_t v2=0; v2 < count; ++v2) {
201                                 buf += string_snprintf("%.5f", hess[v1 * count + v2]);
202                                 if (v1 < count-1 || v2 < count-1) buf += ", ";
203                         }
204                         buf += "\n";
205                 }
206                 buf += ")\n";
207         }
208         if (what & FF_COMPUTE_IHESSIAN) {
209                 buf += string_snprintf("ihess %lux%lu: c(", count, count);
210                 for (size_t v1=0; v1 < count; ++v1) {
211                         for (size_t v2=0; v2 < count; ++v2) {
212                                 buf += string_snprintf("%.5f", ihess[v1 * count + v2]);
213                                 if (v1 < count-1 || v2 < count-1) buf += ", ";
214                         }
215                         buf += "\n";
216                 }
217                 buf += ")\n";
218         }
219         if (what & FF_COMPUTE_HGPROD) {
220                 buf += string_snprintf("ihess %%*%% grad %lu: list(", hgProd.size());
221                 for (size_t px=0; px < hgProd.size(); ++px) {
222                         buf += string_snprintf("c(%d, %d, %d)", hgProd[px].hentry,
223                                                hgProd[px].gentry, hgProd[px].dest);
224                         if (px < hgProd.size() - 1) buf += ", ";
225                 }
226                 buf += ")\n";
227         }
228         mxLogBig(buf);
229 }
230
231 void FitContext::fixHessianSymmetry(int want)
232 {
233         size_t numParam = varGroup->vars.size();
234
235         if (want & (FF_COMPUTE_HESSIAN | FF_COMPUTE_INFO)) {
236                 for (size_t h1=1; h1 < numParam; h1++) {
237                         for (size_t h2=0; h2 < h1; h2++) {
238                                 if (hess[h2 * numParam + h1] != 0) {
239                                         omxRaiseErrorf(globalState, "Hessian/information is not upper triangular");
240                                         break;
241                                 }
242                                 hess[h2 * numParam + h1] = hess[h1 * numParam + h2];
243                         }
244                 }
245         }
246
247         if (want & FF_COMPUTE_IHESSIAN) {
248                 for (size_t h1=1; h1 < numParam; h1++) {
249                         for (size_t h2=0; h2 < h1; h2++) {
250                                 if (ihess[h2 * numParam + h1] != 0) {
251                                         omxRaiseErrorf(globalState, "Inverse Hessian is not upper triangular");
252                                         break;
253                                 }
254                                 ihess[h2 * numParam + h1] = ihess[h1 * numParam + h2];
255                         }
256                 }
257         }
258 }
259
260 static void omxRepopulateRFitFunction(omxFitFunction* oo, double* x, int n)
261 {
262         omxRFitFunction* rFitFunction = (omxRFitFunction*)oo->argStruct;
263
264         SEXP theCall, estimate;
265
266         PROTECT(estimate = allocVector(REALSXP, n));
267         double *est = REAL(estimate);
268         for(int i = 0; i < n ; i++) {
269                 est[i] = x[i];
270         }
271
272         PROTECT(theCall = allocVector(LANGSXP, 4));
273
274         SETCAR(theCall, install("imxUpdateModelValues"));
275         SETCADR(theCall, rFitFunction->model);
276         SETCADDR(theCall, rFitFunction->flatModel);
277         SETCADDDR(theCall, estimate);
278
279         REPROTECT(rFitFunction->model = eval(theCall, R_GlobalEnv), rFitFunction->modelIndex);
280
281         UNPROTECT(2); // theCall, estimate
282 }
283
284 void FitContext::copyParamToModel(omxState* os)
285 {
286         copyParamToModel(os, est);
287 }
288
289 void FitContext::maybeCopyParamToModel(omxState* os)
290 {
291         if (changedEstimates) {
292                 copyParamToModel(os, est);
293                 changedEstimates = false;
294         }
295 }
296
297 void FitContext::copyParamToModel(omxState* os, double *at)
298 {
299         size_t numParam = varGroup->vars.size();
300         if(OMX_DEBUG) {
301                 mxLog("Copying %lu free parameter estimates to model %p", numParam, os);
302         }
303
304         if(numParam == 0) return;
305
306         // Confidence Intervals & Hessian Calculation probe the parameter space
307         // near the best estimate. If stale, we need to restore the best estimate
308         // before returning results to the user.
309         os->stale = at != est;
310
311         os->computeCount++;
312
313         if(OMX_VERBOSE) {
314                 std::string buf;
315                 buf += string_snprintf("Call: %d.%d (%ld) ", os->majorIteration, os->minorIteration, os->computeCount);
316                 buf += ("Estimates: [");
317                 for(size_t k = 0; k < numParam; k++) {
318                         buf += string_snprintf(" %f", at[k]);
319                 }
320                 buf += ("]\n");
321                 mxLogBig(buf);
322         }
323
324         for(size_t k = 0; k < numParam; k++) {
325                 omxFreeVar* freeVar = varGroup->vars[k];
326                 for(size_t l = 0; l < freeVar->locations.size(); l++) {
327                         omxFreeVarLocation *loc = &freeVar->locations[l];
328                         omxMatrix *matrix = os->matrixList[loc->matrix];
329                         int row = loc->row;
330                         int col = loc->col;
331                         omxSetMatrixElement(matrix, row, col, at[k]);
332                         if(OMX_DEBUG) {
333                                 mxLog("Setting location (%d, %d) of matrix %d to value %f for var %lu",
334                                         row, col, loc->matrix, at[k], k);
335                         }
336                 }
337         }
338
339         if (RFitFunction) omxRepopulateRFitFunction(RFitFunction, at, numParam);
340
341         varGroup->markDirty(os);
342
343         if (!os->childList) return;
344
345         for(int i = 0; i < Global->numChildren; i++) {
346                 copyParamToModel(os->childList[i], at);
347         }
348 }
349
350 double *FitContext::take(int want)
351 {
352         if (!(want & (wanted | FF_COMPUTE_ESTIMATE))) {
353                 error("Attempt to take %d but not available", want);
354         }
355
356         double *ret = NULL;
357         switch(want) {
358         case FF_COMPUTE_ESTIMATE:
359                 ret = est;
360                 est = NULL;
361                 break;
362         case FF_COMPUTE_HESSIAN:
363         case FF_COMPUTE_INFO:
364                 ret = hess;
365                 hess = NULL;
366                 break;
367         case FF_COMPUTE_IHESSIAN:
368                 ret = ihess;
369                 ihess = NULL;
370                 break;
371         default:
372                 error("Taking of %d is not implemented", want);
373         }
374         if (!ret) error("Attempt to take %d, already taken", want);
375         return ret;
376 }
377
378 FitContext::~FitContext()
379 {
380         if (est) delete [] est;
381         if (flavor) delete [] flavor;
382         if (grad) delete [] grad;
383         if (hess) delete [] hess;
384         if (ihess) delete [] ihess;
385 }
386
387 omxFitFunction *FitContext::RFitFunction = NULL;
388
389 void FitContext::setRFitFunction(omxFitFunction *rff)
390 {
391         if (rff) {
392                 Global->numThreads = 1;
393                 if (RFitFunction) {
394                         error("You can only create 1 MxRFitFunction per independent model");
395                 }
396         }
397         RFitFunction = rff;
398 }
399
400 Ramsay1975::Ramsay1975(FitContext *fc, int flavor, double caution, int verbose,
401                        double minCaution)
402 {
403         this->fc = fc;
404         this->flavor = flavor;
405         this->verbose = verbose;
406         this->caution = caution;
407         this->minCaution = minCaution;
408         maxCaution = 0.0;
409         highWatermark = std::max(0.5, caution);  // arbitrary guess
410
411         numParam = fc->varGroup->vars.size();
412         prevAdj1.assign(numParam, 0);
413         prevAdj2.resize(numParam);
414         prevEst.resize(numParam);
415         memcpy(prevEst.data(), fc->est, sizeof(double) * numParam);
416 }
417
418 void Ramsay1975::recordEstimate(int px, double newEst)
419 {
420         omxFreeVar *fv = fc->varGroup->vars[px];
421         bool hitBound=false;
422         double param = newEst;
423         if (param < fv->lbound) {
424                 hitBound=true;
425                 param = prevEst[px] - (prevEst[px] - fv->lbound) / 2;
426         }
427         if (param > fv->ubound) {
428                 hitBound=true;
429                 param = prevEst[px] + (fv->ubound - prevEst[px]) / 2;
430         }
431         
432         prevAdj2[px] = prevAdj1[px];
433         prevAdj1[px] = param - prevEst[px];
434         
435         if (verbose >= 4) {
436                 std::string buf;
437                 buf += string_snprintf("~%d~%s: %.4f -> %.4f", px, fv->name, prevEst[px], param);
438                 if (hitBound) {
439                         buf += string_snprintf(" wanted %.4f but hit bound", newEst);
440                 }
441                 if (prevAdj1[px] * prevAdj2[px] < 0) {
442                         buf += " *OSC*";
443                 }
444                 buf += "\n";
445                 mxLogBig(buf);
446         }
447
448         fc->est[px] = param;
449         prevEst[px] = param;
450 }
451
452 void Ramsay1975::apply()
453 {
454         for (size_t px=0; px < numParam; ++px) {
455                 recordEstimate(px, (1 - caution) * fc->est[px] + caution * prevEst[px]);
456         }
457 }
458
459 void Ramsay1975::recalibrate(bool *restart)
460 {
461         double normPrevAdj2 = 0;
462         double normAdjDiff = 0;
463         std::vector<double> adjDiff(numParam);
464
465         // The choice of norm is also arbitrary. Other norms might work better.
466         for (size_t px=0; px < numParam; ++px) {
467                 if (fc->flavor[px] != flavor) continue;
468                 adjDiff[px] = prevAdj1[px] - prevAdj2[px];
469                 normPrevAdj2 += prevAdj2[px] * prevAdj2[px];
470         }
471
472         for (size_t px=0; px < numParam; ++px) {
473                 if (fc->flavor[px] != flavor) continue;
474                 normAdjDiff += adjDiff[px] * adjDiff[px];
475         }
476         if (normAdjDiff == 0) {
477                 return;
478                 //error("Ramsay: no free variables of flavor %d", flavor);
479         }
480
481         double ratio = sqrt(normPrevAdj2 / normAdjDiff);
482         //if (verbose >= 3) mxLog("Ramsay[%d]: sqrt(%.5f/%.5f) = %.5f",
483         // flavor, normPrevAdj2, normAdjDiff, ratio);
484
485         double newCaution = 1 - (1-caution) * ratio;
486         if (newCaution > .95) newCaution = .95;  // arbitrary guess
487         if (newCaution < 0) newCaution /= 2;     // don't get overconfident
488         if (newCaution < minCaution) newCaution = minCaution;
489         if (newCaution < caution) {
490                 caution = newCaution/3 + 2*caution/3;  // don't speed up too fast, arbitrary ratio
491         } else {
492                 caution = newCaution;
493         }
494         maxCaution = std::max(maxCaution, caution);
495         if (caution < highWatermark || (normPrevAdj2 < 1e-3 && normAdjDiff < 1e-3)) {
496                 if (verbose >= 3) mxLog("Ramsay[%d]: %.2f caution", flavor, caution);
497         } else {
498                 if (verbose >= 3) mxLog("Ramsay[%d]: caution %.2f > %.2f, extreme oscillation, restart recommended",
499                                         flavor, caution, highWatermark);
500                 *restart = TRUE;
501         }
502         highWatermark += .02; // arbitrary guess
503 }
504
505 void Ramsay1975::restart()
506 {
507         memcpy(prevEst.data(), fc->est, sizeof(double) * numParam);
508         prevAdj1.assign(numParam, 0);
509         prevAdj2.assign(numParam, 0);
510         highWatermark = 1 - (1 - highWatermark) * .5; // arbitrary guess
511         caution = std::max(caution, highWatermark);   // arbitrary guess
512         maxCaution = std::max(maxCaution, caution);
513         highWatermark = caution;
514         if (verbose >= 3) {
515                 mxLog("Ramsay[%d]: restart with %.2f caution %.2f highWatermark",
516                       flavor, caution, highWatermark);
517         }
518 }
519
520 omxCompute::omxCompute()
521 {
522         varGroup = NULL;
523 }
524
525 omxCompute::~omxCompute()
526 {}
527
528 void omxCompute::initFromFrontend(SEXP rObj)
529 {
530         SEXP slotValue;
531         PROTECT(slotValue = GET_SLOT(rObj, install("id")));
532         if (length(slotValue) == 1) {
533                 int id = INTEGER(slotValue)[0];
534                 varGroup = Global->findVarGroup(id);
535         }
536
537         if (!varGroup) {
538                 if (!R_has_slot(rObj, install("free.set"))) {
539                         varGroup = Global->freeGroup[0];
540                 } else {
541                         PROTECT(slotValue = GET_SLOT(rObj, install("free.set")));
542                         if (length(slotValue) != 0) {
543                                 // it's a free.set with no free variables
544                                 varGroup = Global->findVarGroup(-1);
545                         } else {
546                                 varGroup = Global->freeGroup[0];
547                         }
548                 }
549         }
550 }
551
552 class omxComputeSequence : public omxCompute {
553         typedef omxCompute super;
554         std::vector< omxCompute* > clist;
555
556  public:
557         virtual void initFromFrontend(SEXP rObj);
558         virtual void compute(FitContext *fc);
559         virtual void reportResults(FitContext *fc, MxRList *out);
560         virtual double getOptimizerStatus();
561         virtual ~omxComputeSequence();
562 };
563
564 class omxComputeIterate : public omxCompute {
565         typedef omxCompute super;
566         std::vector< omxCompute* > clist;
567         int maxIter;
568         double tolerance;
569         int verbose;
570
571  public:
572         virtual void initFromFrontend(SEXP rObj);
573         virtual void compute(FitContext *fc);
574         virtual void reportResults(FitContext *fc, MxRList *out);
575         virtual double getOptimizerStatus();
576         virtual ~omxComputeIterate();
577 };
578
579 class omxComputeOnce : public omxCompute {
580         typedef omxCompute super;
581         std::vector< omxMatrix* > algebras;
582         std::vector< omxExpectation* > expectations;
583         int verbose;
584         const char *context;
585         bool mac;
586         bool fit;
587         bool gradient;
588         bool hessian;
589         bool ihessian;
590         bool infoMat;
591         bool hgprod;
592
593  public:
594         virtual void initFromFrontend(SEXP rObj);
595         virtual omxFitFunction *getFitFunction();
596         virtual void compute(FitContext *fc);
597         virtual void reportResults(FitContext *fc, MxRList *out);
598 };
599
600 class ComputeEM : public omxCompute {
601         typedef omxCompute super;
602         std::vector< omxExpectation* > expectations;
603         omxCompute *fit1;
604         omxCompute *fit2;
605         int maxIter;
606         int mstepIter;
607         int totalMstepIter;
608         double tolerance;
609         int verbose;
610         bool useRamsay;
611         bool information;
612         std::vector<Ramsay1975*> ramsay;
613         std::vector<double*> estHistory;
614         FitContext *recentFC;  //nice if can use std::unique_ptr
615         std::vector<double> optimum;
616         std::vector<double> stdError;
617         double bestFit;
618         static const double MIDDLE_START = 0.21072103131565256273; // -log(.9)*2 constexpr
619         static const double MIDDLE_END = 0.0020010006671670687271; // -log(.999)*2
620
621         void setExpectationContext(const char *context);
622
623  public:
624         virtual void initFromFrontend(SEXP rObj);
625         virtual void compute(FitContext *fc);
626         virtual void reportResults(FitContext *fc, MxRList *out);
627         virtual double getOptimizerStatus();
628         virtual ~ComputeEM();
629 };
630
631 static class omxCompute *newComputeSequence()
632 { return new omxComputeSequence(); }
633
634 static class omxCompute *newComputeIterate()
635 { return new omxComputeIterate(); }
636
637 static class omxCompute *newComputeOnce()
638 { return new omxComputeOnce(); }
639
640 static class omxCompute *newComputeEM()
641 { return new ComputeEM(); }
642
643 struct omxComputeTableEntry {
644         char name[32];
645         omxCompute *(*ctor)();
646 };
647
648 static const struct omxComputeTableEntry omxComputeTable[] = {
649         {"MxComputeEstimatedHessian", &newComputeEstimatedHessian},
650         {"MxComputeGradientDescent", &newComputeGradientDescent},
651         {"MxComputeSequence", &newComputeSequence },
652         {"MxComputeIterate", &newComputeIterate },
653         {"MxComputeOnce", &newComputeOnce },
654         {"MxComputeNewtonRaphson", &newComputeNewtonRaphson},
655         {"MxComputeEM", &newComputeEM },
656 };
657
658 omxCompute *omxNewCompute(omxState* os, const char *type)
659 {
660         omxCompute *got = NULL;
661
662         for (size_t fx=0; fx < OMX_STATIC_ARRAY_SIZE(omxComputeTable); fx++) {
663                 const struct omxComputeTableEntry *entry = omxComputeTable + fx;
664                 if(strcmp(type, entry->name) == 0) {
665                         got = entry->ctor();
666                         break;
667                 }
668         }
669
670         if (!got) error("Compute %s is not implemented", type);
671
672         return got;
673 }
674
675 void omxComputeSequence::initFromFrontend(SEXP rObj)
676 {
677         super::initFromFrontend(rObj);
678
679         SEXP slotValue;
680         PROTECT(slotValue = GET_SLOT(rObj, install("steps")));
681
682         for (int cx = 0; cx < length(slotValue); cx++) {
683                 SEXP step = VECTOR_ELT(slotValue, cx);
684                 SEXP s4class;
685                 PROTECT(s4class = STRING_ELT(getAttrib(step, install("class")), 0));
686                 omxCompute *compute = omxNewCompute(globalState, CHAR(s4class));
687                 compute->initFromFrontend(step);
688                 if (isErrorRaised(globalState)) break;
689                 clist.push_back(compute);
690         }
691 }
692
693 void omxComputeSequence::compute(FitContext *fc)
694 {
695         for (size_t cx=0; cx < clist.size(); ++cx) {
696                 FitContext *context = fc;
697                 if (fc->varGroup != clist[cx]->varGroup) {
698                         context = new FitContext(fc, clist[cx]->varGroup);
699                 }
700                 clist[cx]->compute(context);
701                 if (context != fc) context->updateParentAndFree();
702                 if (isErrorRaised(globalState)) break;
703         }
704 }
705
706 void omxComputeSequence::reportResults(FitContext *fc, MxRList *out)
707 {
708         // put this stuff in a new list?
709         // merge with Iterate TODO
710         for (size_t cx=0; cx < clist.size(); ++cx) {
711                 FitContext *context = fc;
712                 if (fc->varGroup != clist[cx]->varGroup) {
713                         context = new FitContext(fc, clist[cx]->varGroup);
714                 }
715                 clist[cx]->reportResults(context, out);
716                 if (context != fc) context->updateParentAndFree();
717                 if (isErrorRaised(globalState)) break;
718         }
719 }
720
721 double omxComputeSequence::getOptimizerStatus()
722 {
723         // for backward compatibility, not indended to work generally
724         for (size_t cx=0; cx < clist.size(); ++cx) {
725                 double got = clist[cx]->getOptimizerStatus();
726                 if (got != NA_REAL) return got;
727         }
728         return NA_REAL;
729 }
730
731 omxComputeSequence::~omxComputeSequence()
732 {
733         for (size_t cx=0; cx < clist.size(); ++cx) {
734                 delete clist[cx];
735         }
736 }
737
738 void omxComputeIterate::initFromFrontend(SEXP rObj)
739 {
740         SEXP slotValue;
741
742         super::initFromFrontend(rObj);
743
744         PROTECT(slotValue = GET_SLOT(rObj, install("maxIter")));
745         maxIter = INTEGER(slotValue)[0];
746
747         PROTECT(slotValue = GET_SLOT(rObj, install("tolerance")));
748         tolerance = REAL(slotValue)[0];
749         if (tolerance <= 0) error("tolerance must be positive");
750
751         PROTECT(slotValue = GET_SLOT(rObj, install("steps")));
752
753         for (int cx = 0; cx < length(slotValue); cx++) {
754                 SEXP step = VECTOR_ELT(slotValue, cx);
755                 SEXP s4class;
756                 PROTECT(s4class = STRING_ELT(getAttrib(step, install("class")), 0));
757                 omxCompute *compute = omxNewCompute(globalState, CHAR(s4class));
758                 compute->initFromFrontend(step);
759                 if (isErrorRaised(globalState)) break;
760                 clist.push_back(compute);
761         }
762
763         PROTECT(slotValue = GET_SLOT(rObj, install("verbose")));
764         verbose = asInteger(slotValue);
765 }
766
767 void omxComputeIterate::compute(FitContext *fc)
768 {
769         int iter = 0;
770         double prevFit = 0;
771         double mac = tolerance * 10;
772         while (1) {
773                 for (size_t cx=0; cx < clist.size(); ++cx) {
774                         FitContext *context = fc;
775                         if (fc->varGroup != clist[cx]->varGroup) {
776                                 context = new FitContext(fc, clist[cx]->varGroup);
777                         }
778                         clist[cx]->compute(context);
779                         if (context != fc) context->updateParentAndFree();
780                         if (isErrorRaised(globalState)) break;
781                 }
782                 if (fc->wanted & FF_COMPUTE_MAXABSCHANGE) {
783                         if (fc->mac < 0) {
784                                 warning("MAC estimated at %.4f; something is wrong", fc->mac);
785                                 break;
786                         } else {
787                                 mac = fc->mac;
788                                 if (verbose) mxLog("ComputeIterate: mac %.9g", mac);
789                         }
790                 }
791                 if (fc->wanted & FF_COMPUTE_FIT) {
792                         if (fc->fit == 0) {
793                                 warning("Fit estimated at 0; something is wrong");
794                                 break;
795                         }
796                         if (prevFit != 0) {
797                                 double change = prevFit - fc->fit;
798                                 if (verbose) mxLog("ComputeIterate: fit %.9g change %.9g", fc->fit, change);
799                                 mac = fabs(change);
800                         } else {
801                                 if (verbose) mxLog("ComputeIterate: initial fit %.9g", fc->fit);
802                         }
803                         prevFit = fc->fit;
804                 }
805                 if (!(fc->wanted & (FF_COMPUTE_MAXABSCHANGE | FF_COMPUTE_FIT))) {
806                         omxRaiseErrorf(globalState, "ComputeIterate: neither MAC nor fit available");
807                 }
808                 if (isErrorRaised(globalState) || ++iter > maxIter || mac < tolerance) break;
809         }
810 }
811
812 void omxComputeIterate::reportResults(FitContext *fc, MxRList *out)
813 {
814         for (size_t cx=0; cx < clist.size(); ++cx) {
815                 FitContext *context = fc;
816                 if (fc->varGroup != clist[cx]->varGroup) {
817                         context = new FitContext(fc, clist[cx]->varGroup);
818                 }
819                 clist[cx]->reportResults(context, out);
820                 if (context != fc) context->updateParentAndFree();
821                 if (isErrorRaised(globalState)) break;
822         }
823 }
824
825 double omxComputeIterate::getOptimizerStatus()
826 {
827         // for backward compatibility, not indended to work generally
828         for (size_t cx=0; cx < clist.size(); ++cx) {
829                 double got = clist[cx]->getOptimizerStatus();
830                 if (got != NA_REAL) return got;
831         }
832         return NA_REAL;
833 }
834
835 omxComputeIterate::~omxComputeIterate()
836 {
837         for (size_t cx=0; cx < clist.size(); ++cx) {
838                 delete clist[cx];
839         }
840 }
841
842 void ComputeEM::initFromFrontend(SEXP rObj)
843 {
844         recentFC = NULL;
845
846         SEXP slotValue;
847         SEXP s4class;
848
849         super::initFromFrontend(rObj);
850
851         PROTECT(slotValue = GET_SLOT(rObj, install("maxIter")));
852         maxIter = INTEGER(slotValue)[0];
853
854         PROTECT(slotValue = GET_SLOT(rObj, install("information")));
855         information = asLogical(slotValue);
856
857         PROTECT(slotValue = GET_SLOT(rObj, install("ramsay")));
858         useRamsay = asLogical(slotValue);
859
860         PROTECT(slotValue = GET_SLOT(rObj, install("tolerance")));
861         tolerance = REAL(slotValue)[0];
862         if (tolerance <= 0) error("tolerance must be positive");
863
864         PROTECT(slotValue = GET_SLOT(rObj, install("what")));
865         for (int wx=0; wx < length(slotValue); ++wx) {
866                 int objNum = INTEGER(slotValue)[wx];
867                 omxExpectation *expectation = globalState->expectationList[objNum];
868                 setFreeVarGroup(expectation, varGroup);
869                 omxCompleteExpectation(expectation);
870                 expectations.push_back(expectation);
871         }
872
873         PROTECT(slotValue = GET_SLOT(rObj, install("mstep.fit")));
874         PROTECT(s4class = STRING_ELT(getAttrib(slotValue, install("class")), 0));
875         fit1 = omxNewCompute(globalState, CHAR(s4class));
876         fit1->initFromFrontend(slotValue);
877
878         PROTECT(slotValue = GET_SLOT(rObj, install("fit")));
879         PROTECT(s4class = STRING_ELT(getAttrib(slotValue, install("class")), 0));
880         fit2 = omxNewCompute(globalState, CHAR(s4class));
881         fit2->initFromFrontend(slotValue);
882
883         PROTECT(slotValue = GET_SLOT(rObj, install("verbose")));
884         verbose = asInteger(slotValue);
885 }
886
887 void ComputeEM::setExpectationContext(const char *context)
888 {
889         for (size_t wx=0; wx < expectations.size(); ++wx) {
890                 omxExpectation *expectation = expectations[wx];
891                 if (verbose >= 4) mxLog("ComputeEM: expectation[%lu] %s context %s", wx, expectation->name, context);
892                 omxExpectationCompute(expectation, context);
893         }
894 }
895
896 void pda(const double *ar, int rows, int cols);
897
898 void ComputeEM::compute(FitContext *fc)
899 {
900         int totalMstepIter = 0;
901         int iter = 0;
902         double prevFit = 0;
903         double mac = tolerance * 10;
904         bool converged = false;
905         const size_t freeVars = fc->varGroup->vars.size();
906         const int freeVarsEM = (int) fit1->varGroup->vars.size();
907         bool in_middle = false;
908
909         OMXZERO(fc->flavor, freeVars);
910
911         FitContext *tmp = new FitContext(fc, fit1->varGroup);
912         for (int vx=0; vx < freeVarsEM; ++vx) {
913                 fc->flavor[tmp->mapToParent[vx]] = 1;
914         }
915         tmp->updateParentAndFree();
916
917         ramsay.push_back(new Ramsay1975(fc, int(ramsay.size()), 0, verbose, -1.25)); // other param
918         ramsay.push_back(new Ramsay1975(fc, int(ramsay.size()), 0, verbose, -1));    // EM param
919
920         if (verbose >= 1) mxLog("ComputeEM: Welcome, tolerance=%g ramsay=%d info=%d flavors=%ld",
921                                 tolerance, useRamsay, information, ramsay.size());
922
923         while (1) {
924                 setExpectationContext("EM");
925
926                 if (recentFC) delete recentFC;
927                 recentFC = new FitContext(fc, fit1->varGroup);
928                 fit1->compute(recentFC);
929                 mstepIter = recentFC->iterations;
930                 recentFC->updateParent();
931
932                 setExpectationContext("");
933
934                 {
935                         FitContext *context = fc;
936                         if (fc->varGroup != fit2->varGroup) {
937                                 context = new FitContext(fc, fit2->varGroup);
938                         }
939
940                         // For IFA, PREOPTIMIZE updates latent distribution parameters
941                         omxFitFunction *ff2 = fit2->getFitFunction();
942                         if (ff2) omxFitFunctionCompute(ff2, FF_COMPUTE_PREOPTIMIZE, context);
943
944                         if (!useRamsay) {
945                                 fc->maybeCopyParamToModel(globalState);
946                         } else {
947                                 context->updateParent();
948
949                                 bool wantRestart;
950                                 if (iter > 3 && iter % 3 == 0) {
951                                         for (size_t rx=0; rx < ramsay.size(); ++rx) {
952                                                 ramsay[rx]->recalibrate(&wantRestart);
953                                         }
954                                 }
955                                 for (size_t rx=0; rx < ramsay.size(); ++rx) {
956                                         ramsay[rx]->apply();
957                                 }
958                                 fc->copyParamToModel(globalState);
959                         }
960
961                         fit2->compute(context);
962                         if (context != fc) context->updateParentAndFree();
963                 }
964
965                 totalMstepIter += mstepIter;
966
967                 if (!(fc->wanted & FF_COMPUTE_FIT)) {
968                         omxRaiseErrorf(globalState, "ComputeEM: fit not available");
969                         break;
970                 }
971                 if (fc->fit == 0) {
972                         omxRaiseErrorf(globalState, "Fit estimated at 0; something is wrong");
973                         break;
974                 }
975                 double change = 0;
976                 if (prevFit != 0) {
977                         change = prevFit - fc->fit;
978                         if (0 < change && change < MIDDLE_START) in_middle = true;
979                         if (verbose >= 2) mxLog("ComputeEM[%d]: iter %d fit %.9g change %.9g",
980                                                 iter, mstepIter, fc->fit, change);
981                         mac = fabs(change);
982                 } else {
983                         if (verbose >= 2) mxLog("ComputeEM: iter %d initial fit %.9g",
984                                                 mstepIter, fc->fit);
985                 }
986
987                 if (in_middle && change > MIDDLE_END) estHistory.push_back(recentFC->take(FF_COMPUTE_ESTIMATE));
988                 prevFit = fc->fit;
989                 converged = mac < tolerance;
990                 if (isErrorRaised(globalState) || ++iter > maxIter || converged) break;
991         }
992
993         bestFit = fc->fit;
994         if (verbose >= 1) mxLog("ComputeEM: cycles %d/%d total mstep %d fit %f",
995                                 iter, maxIter,totalMstepIter, bestFit);
996
997         if (!converged || !information) return;
998
999         if (estHistory.size() < 2) {
1000                 if (verbose >= 1) mxLog("ComputeEM: history too short to estimate SEs; try increasing EM tolerance");
1001                 return;
1002         }
1003
1004         // what about latent distribution parameters? TODO
1005
1006         recentFC->fixHessianSymmetry(FF_COMPUTE_IHESSIAN);
1007         double *ihess = recentFC->take(FF_COMPUTE_IHESSIAN);
1008
1009         double semTolerance = sqrt(tolerance);
1010         optimum.resize(freeVars);
1011         memcpy(optimum.data(), fc->est, sizeof(double) * freeVars);
1012         omxBuffer<double> rij(freeVarsEM * freeVarsEM);
1013         setExpectationContext("EM");
1014
1015         // could parallelize by variable instead of at lower level, not sure which is better TODO
1016         for (int vx=0; vx < freeVarsEM && !isErrorRaised(globalState); ++vx) {
1017                 int base = vx * freeVarsEM;
1018                 bool semConverged;
1019                 for (size_t cx = 0; cx < estHistory.size(); ++cx) {
1020                         if (verbose >= 2) mxLog("ComputeEM: probing param %d from history %ld/%ld",
1021                                                 vx, cx, estHistory.size());
1022
1023                         memcpy(fc->est, optimum.data(), sizeof(double) * freeVars);
1024                         FitContext *emfc = new FitContext(fc, fit1->varGroup);
1025                         emfc->est[vx] = estHistory[cx][vx];
1026                         emfc->copyParamToModel(globalState);
1027                         fit1->compute(emfc);
1028
1029                         double denom = estHistory[cx][vx] - optimum[emfc->mapToParent[vx]];
1030                         if (verbose >= 1 && fabs(denom) < 1e-5) {
1031                                 mxLog("ComputeEM: param %d history %ld denom=%f < 1e-5 (this is bad)", vx, cx, denom);
1032                         }
1033
1034                         semConverged = true;
1035                         for (int v1=0; v1 < freeVarsEM; ++v1) {
1036                                 double got = (emfc->est[v1] - optimum[emfc->mapToParent[v1]]) / denom;
1037                                 if (fabs(got - rij[base + v1]) >= semTolerance) semConverged = false;
1038                                 rij[base + v1] = got;
1039                         }
1040                         //pda(rij.data() + base, 1, freeVarsEM);
1041                         if (verbose >= 3) mxLog("rij[%d]=%f", vx, rij[base + vx]); // useful?
1042                         delete emfc;
1043                         if ((cx > 0 && semConverged) || isErrorRaised(globalState)) break;
1044                 }
1045         }
1046
1047         memcpy(fc->est, optimum.data(), sizeof(double) * freeVars);
1048         fc->copyParamToModel(globalState);
1049
1050         //pda(rij.data(), freeVarsEM, freeVarsEM);
1051
1052         // rij = I-rij
1053         for (int v1=0; v1 < freeVarsEM; ++v1) {
1054                 for (int v2=0; v2 < freeVarsEM; ++v2) {
1055                         int cell = v1 * freeVarsEM + v2;
1056                         double entry = rij[cell];
1057                         if (v1 == v2) entry = 1 - entry;
1058                         else entry = -entry;
1059                         rij[cell] = entry;
1060                 }
1061         }
1062         // make symmetric
1063         for (int v1=1; v1 < freeVarsEM; ++v1) {
1064                 for (int v2=0; v2 < v1; ++v2) {
1065                         int c1 = v1 * freeVarsEM + v2;
1066                         int c2 = v2 * freeVarsEM + v1;
1067                         double mean = (rij[c1] + rij[c2])/2;
1068                         rij[c1] = mean;
1069                         rij[c2] = mean;
1070                 }
1071         }
1072
1073         //mxLog("symm");
1074         //pda(rij.data(), freeVarsEM, freeVarsEM);
1075
1076         //pda(ihess, freeVarsEM, freeVarsEM);
1077
1078         // ihess = ihess %*% rij^{-1}
1079         omxBuffer<int> ipiv(freeVarsEM);
1080         int info;
1081         F77_CALL(dgesv)(&freeVarsEM, &freeVarsEM, rij.data(), &freeVarsEM,  // dsysv? TODO
1082                         ipiv.data(), ihess, &freeVarsEM, &info);
1083         if (info < 0) error("dgesv %d", info);
1084         if (info > 0) {
1085                 omxRaiseErrorf(globalState, "EM map is not positive definite %d", info);
1086                 return;
1087         }
1088
1089         //pda(ihess, freeVarsEM, freeVarsEM);
1090
1091         stdError.resize(freeVars);
1092         for (int v1=0; v1 < freeVarsEM; ++v1) {
1093                 int cell = v1 * freeVarsEM + v1;
1094                 stdError[ recentFC->mapToParent[v1] ] = sqrt(ihess[cell]);
1095         }       
1096
1097         delete [] ihess;
1098 }
1099
1100 void ComputeEM::reportResults(FitContext *fc, MxRList *out)
1101 {
1102         out->push_back(std::make_pair(mkChar("minimum"), ScalarReal(bestFit)));
1103         out->push_back(std::make_pair(mkChar("Minus2LogLikelihood"), ScalarReal(bestFit)));
1104
1105         size_t numFree = fc->varGroup->vars.size();
1106         if (!numFree) return;
1107
1108         if (optimum.size() == numFree) {
1109                 SEXP Rvec;
1110                 PROTECT(Rvec = allocVector(REALSXP, numFree));
1111                 memcpy(REAL(Rvec), optimum.data(), sizeof(double)*numFree);
1112                 out->push_back(std::make_pair(mkChar("estimate"), Rvec));
1113         }
1114
1115         if (stdError.size() == numFree) {
1116                 // make conditional TODO
1117                 SEXP Rvec;
1118                 PROTECT(Rvec = allocVector(REALSXP, numFree));
1119                 memcpy(REAL(Rvec), stdError.data(), sizeof(double)*numFree);
1120                 out->push_back(std::make_pair(mkChar("standardErrors"), Rvec));
1121         }
1122 }
1123
1124 double ComputeEM::getOptimizerStatus()
1125 {
1126         // for backward compatibility, not indended to work generally
1127         return NA_REAL;
1128 }
1129
1130 ComputeEM::~ComputeEM()
1131 {
1132         for (size_t rx=0; rx < ramsay.size(); ++rx) {
1133                 delete ramsay[rx];
1134         }
1135         ramsay.clear();
1136
1137         delete fit1;
1138         delete fit2;
1139
1140         for (size_t hx=0; hx < estHistory.size(); ++hx) {
1141                 delete [] estHistory[hx];
1142         }
1143         estHistory.clear();
1144         if (recentFC) delete recentFC;
1145 }
1146
1147 void omxComputeOnce::initFromFrontend(SEXP rObj)
1148 {
1149         super::initFromFrontend(rObj);
1150
1151         SEXP slotValue;
1152         PROTECT(slotValue = GET_SLOT(rObj, install("what")));
1153         for (int wx=0; wx < length(slotValue); ++wx) {
1154                 int objNum = INTEGER(slotValue)[wx];
1155                 if (objNum >= 0) {
1156                         omxMatrix *algebra = globalState->algebraList[objNum];
1157                         if (algebra->fitFunction) {
1158                                 setFreeVarGroup(algebra->fitFunction, varGroup);
1159                                 omxCompleteFitFunction(algebra);
1160                         }
1161                         algebras.push_back(algebra);
1162                 } else {
1163                         omxExpectation *expectation = globalState->expectationList[~objNum];
1164                         setFreeVarGroup(expectation, varGroup);
1165                         omxCompleteExpectation(expectation);
1166                         expectations.push_back(expectation);
1167                 }
1168         }
1169
1170         PROTECT(slotValue = GET_SLOT(rObj, install("verbose")));
1171         verbose = asInteger(slotValue);
1172
1173         context = "";
1174
1175         PROTECT(slotValue = GET_SLOT(rObj, install("context")));
1176         if (length(slotValue) == 0) {
1177                 // OK
1178         } else if (length(slotValue) == 1) {
1179                 SEXP elem;
1180                 PROTECT(elem = STRING_ELT(slotValue, 0));
1181                 context = CHAR(elem);
1182         }
1183
1184         PROTECT(slotValue = GET_SLOT(rObj, install("maxAbsChange")));
1185         mac = asLogical(slotValue);
1186
1187         PROTECT(slotValue = GET_SLOT(rObj, install("fit")));
1188         fit = asLogical(slotValue);
1189
1190         PROTECT(slotValue = GET_SLOT(rObj, install("gradient")));
1191         gradient = asLogical(slotValue);
1192
1193         PROTECT(slotValue = GET_SLOT(rObj, install("hessian")));
1194         hessian = asLogical(slotValue);
1195
1196         PROTECT(slotValue = GET_SLOT(rObj, install("information")));
1197         infoMat = asLogical(slotValue);
1198
1199         if (hessian && infoMat) error("Cannot compute the Hessian and Fisher Information matrix simultaneously");
1200
1201         PROTECT(slotValue = GET_SLOT(rObj, install("ihessian")));
1202         ihessian = asLogical(slotValue);
1203
1204         PROTECT(slotValue = GET_SLOT(rObj, install("hgprod")));
1205         hgprod = asLogical(slotValue);
1206
1207         if (algebras.size() == 1 && algebras[0]->fitFunction) {
1208                 omxFitFunction *ff = algebras[0]->fitFunction;
1209                 if (gradient && !ff->gradientAvailable) {
1210                         error("Gradient requested but not available");
1211                 }
1212                 if ((hessian || ihessian || hgprod) && !ff->hessianAvailable) {
1213                         // add a separate flag for hgprod TODO
1214                         error("Hessian requested but not available");
1215                 }
1216                 // add check for information TODO
1217         }
1218 }
1219
1220 omxFitFunction *omxComputeOnce::getFitFunction()
1221 {
1222         if (algebras.size() == 1 && algebras[0]->fitFunction) {
1223                 return algebras[0]->fitFunction;
1224         } else {
1225                 return NULL;
1226         }
1227 }
1228
1229 void omxComputeOnce::compute(FitContext *fc)
1230 {
1231         if (algebras.size()) {
1232                 int want = 0;
1233                 size_t numParam = fc->varGroup->vars.size();
1234                 if (mac) {
1235                         want |= FF_COMPUTE_MAXABSCHANGE;
1236                         fc->mac = 0;
1237                 }
1238                 if (fit) {
1239                         want |= FF_COMPUTE_FIT;
1240                         fc->fit = 0;
1241                 }
1242                 if (gradient) {
1243                         want |= FF_COMPUTE_GRADIENT;
1244                         OMXZERO(fc->grad, numParam);
1245                 }
1246                 if (hessian) {
1247                         want |= FF_COMPUTE_HESSIAN;
1248                         OMXZERO(fc->hess, numParam * numParam);
1249                 }
1250                 if (infoMat) {
1251                         want |= FF_COMPUTE_INFO;
1252                         OMXZERO(fc->hess, numParam * numParam);
1253                 }
1254                 if (ihessian) {
1255                         want |= FF_COMPUTE_IHESSIAN;
1256                         OMXZERO(fc->ihess, numParam * numParam);
1257                 }
1258                 if (hgprod) {
1259                         want |= FF_COMPUTE_HGPROD;
1260                         fc->hgProd.resize(0);
1261                 }
1262                 if (!want) return;
1263
1264                 for (size_t wx=0; wx < algebras.size(); ++wx) {
1265                         omxMatrix *algebra = algebras[wx];
1266                         if (algebra->fitFunction) {
1267                                 if (verbose) mxLog("ComputeOnce: fit %p want %d",
1268                                                    algebra->fitFunction, want);
1269
1270                                 omxFitFunctionCompute(algebra->fitFunction, FF_COMPUTE_PREOPTIMIZE, fc);
1271                                 fc->maybeCopyParamToModel(globalState);
1272
1273                                 omxFitFunctionCompute(algebra->fitFunction, want, fc);
1274                                 fc->fit = algebra->data[0];
1275                                 fc->fixHessianSymmetry(want);
1276                         } else {
1277                                 if (verbose) mxLog("ComputeOnce: algebra %p", algebra);
1278                                 omxForceCompute(algebra);
1279                         }
1280                 }
1281         } else if (expectations.size()) {
1282                 for (size_t wx=0; wx < expectations.size(); ++wx) {
1283                         omxExpectation *expectation = expectations[wx];
1284                         if (verbose) mxLog("ComputeOnce: expectation[%lu] %p context %s", wx, expectation, context);
1285                         omxExpectationCompute(expectation, context);
1286                 }
1287         }
1288 }
1289
1290 void omxComputeOnce::reportResults(FitContext *fc, MxRList *out)
1291 {
1292         if (algebras.size()==0 || algebras[0]->fitFunction == NULL) return;
1293
1294         omxMatrix *algebra = algebras[0];
1295
1296         omxPopulateFitFunction(algebra, out);
1297
1298         out->push_back(std::make_pair(mkChar("minimum"), ScalarReal(fc->fit)));
1299         out->push_back(std::make_pair(mkChar("Minus2LogLikelihood"), ScalarReal(fc->fit)));
1300
1301         size_t numFree = fc->varGroup->vars.size();
1302         if (numFree) {
1303                 SEXP estimate;
1304                 PROTECT(estimate = allocVector(REALSXP, numFree));
1305                 memcpy(REAL(estimate), fc->est, sizeof(double)*numFree);
1306                 out->push_back(std::make_pair(mkChar("estimate"), estimate));
1307
1308                 if (gradient) {
1309                         SEXP Rgradient;
1310                         PROTECT(Rgradient = allocVector(REALSXP, numFree));
1311                         memcpy(REAL(Rgradient), fc->grad, sizeof(double) * numFree);
1312                         out->push_back(std::make_pair(mkChar("gradient"), Rgradient));
1313                 }
1314
1315                 if (hessian) {
1316                         SEXP Rhessian;
1317                         PROTECT(Rhessian = allocMatrix(REALSXP, numFree, numFree));
1318                         memcpy(REAL(Rhessian), fc->hess, sizeof(double) * numFree * numFree);
1319                         out->push_back(std::make_pair(mkChar("hessian"), Rhessian));
1320                 }
1321
1322                 if (infoMat) {
1323                         SEXP Rhessian;
1324                         PROTECT(Rhessian = allocMatrix(REALSXP, numFree, numFree));
1325                         memcpy(REAL(Rhessian), fc->hess, sizeof(double) * numFree * numFree);
1326                         out->push_back(std::make_pair(mkChar("information"), Rhessian));
1327                 }
1328
1329                 if (ihessian) {
1330                         SEXP Rihessian;
1331                         PROTECT(Rihessian = allocMatrix(REALSXP, numFree, numFree));
1332                         memcpy(REAL(Rihessian), fc->ihess, sizeof(double) * numFree * numFree);
1333                         out->push_back(std::make_pair(mkChar("ihessian"), Rihessian));
1334                 }
1335
1336                 if (hgprod) {
1337                         // TODO
1338                 }
1339         }
1340 }