Ignore non-convergence in ComputeEM M-step
[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 <algorithm>
18 #include <stdarg.h>
19
20 #include "Compute.h"
21 #include "Eigen/Cholesky"
22 #include "omxState.h"
23 #include "omxExportBackendState.h"
24 #include "omxRFitFunction.h"
25 #include "matrix.h"
26 #include "omxBuffer.h"
27 #include "omxState.h"
28
29 void pda(const double *ar, int rows, int cols);
30
31 // static bool compareBlocks(const HessianBlock*& lhs, const HessianBlock*& rhs)
32 // {
33 //      return lhs->vars[0] < rhs->vars[0];
34 // }
35
36 void FitContext::queue(HessianBlock *hb)
37 {
38         if (hb->vars.size() == 0) {
39                 delete hb;
40                 return;
41         }
42
43         allBlocks.push_back(hb);
44         //std::push_heap(allBlocks.begin(), allBlocks.end(), compareBlocks); // maybe TODO
45 }
46
47 void FitContext::negateHessian()
48 {
49         // this assumes that we haven't processed the blocks yet, need to check TODO
50         for (size_t bx=0; bx < allBlocks.size(); ++bx) {
51                 allBlocks[bx]->mat *= -1.0;
52         }
53 }
54
55 void FitContext::refreshDenseHess()
56 {
57         if (haveDenseHess) return;
58
59         hess.triangularView<Eigen::Upper>().setZero();
60
61         for (size_t bx=0; bx < allBlocks.size(); ++bx) {
62                 HessianBlock *hb = allBlocks[bx];
63
64                 std::vector<int> &map = hb->vars;
65                 size_t bsize = map.size();
66
67                 for (size_t v1=0; v1 < bsize; ++v1) {
68                         for (size_t v2=0; v2 <= v1; ++v2) {
69                                 hess(map[v2], map[v1]) += hb->mat(v2, v1);
70                         }
71                 }
72         }
73
74         haveDenseHess = true;
75 }
76
77 void FitContext::copyDenseHess(double *dest)
78 {
79         refreshDenseHess();
80         for (size_t v1=0; v1 < numParam; ++v1) {
81                 for (size_t v2=0; v2 <= v1; ++v2) {
82                         double coef = hess.selfadjointView<Eigen::Upper>()(v2,v1);
83                         if (v1==v2) {
84                                 dest[v1 * numParam + v2] = coef;
85                         } else {
86                                 dest[v1 * numParam + v2] = coef;
87                                 dest[v2 * numParam + v1] = coef;
88                         }
89                 }
90         }
91 }
92
93 double *FitContext::getDenseHessUninitialized()
94 {
95         // Assume the caller is going to fill it out
96         haveDenseHess = true;
97         haveDenseIHess = false;
98         return hess.data();
99 }
100
101 void FitContext::refreshDenseIHess()
102 {
103         if (haveDenseIHess) return;
104
105         refreshDenseHess();
106         ihess = hess;
107         Matrix wmat(ihess.data(), numParam, numParam);
108         InvertSymmetricIndef(wmat, 'U');
109
110         haveDenseIHess = true;
111 }
112
113 Eigen::VectorXd FitContext::ihessGradProd()
114 {
115         refreshDenseIHess();
116         return ihess.selfadjointView<Eigen::Upper>() * grad;
117 }
118
119 Eigen::VectorXd FitContext::ihessDiag()
120 {
121         refreshDenseIHess();
122         return ihess.diagonal();
123 }
124
125 double *FitContext::getDenseIHessUninitialized()
126 {
127         // Assume the caller is going to fill it out
128         haveDenseIHess = true;
129         haveDenseHess = false;
130         return ihess.data();
131 }
132
133 Eigen::MatrixXd &FitContext::getDenseIHess()
134 {
135         refreshDenseIHess();
136         return ihess;
137 }
138
139 void FitContext::copyDenseIHess(double *dest)
140 {
141         refreshDenseIHess();
142         for (size_t v1=0; v1 < numParam; ++v1) {
143                 for (size_t v2=0; v2 <= v1; ++v2) {
144                         double coef = ihess.selfadjointView<Eigen::Upper>()(v2,v1);
145                         if (v1==v2) {
146                                 dest[v1 * numParam + v2] = coef;
147                         } else {
148                                 dest[v1 * numParam + v2] = coef;
149                                 dest[v2 * numParam + v1] = coef;
150                         }
151                 }
152         }
153 }
154
155 double *FitContext::getDenseHessianish()
156 {
157         if (haveDenseHess) return hess.data();
158         if (haveDenseIHess) return ihess.data();
159         // try harder TODO
160         return NULL;
161 }
162
163 HessianBlock *HessianBlock::clone()
164 {
165         HessianBlock *hb = new HessianBlock;
166         hb->vars = vars;
167         hb->mat.resize(vars.size(), vars.size());
168         return hb;
169 }
170
171 bool HessianBlock::posDefinite()
172 {
173         Eigen::LLT<Eigen::MatrixXd> llt;
174         llt.compute(mat);
175         return llt.info() == Eigen::Success;
176 }
177
178 void FitContext::init()
179 {
180         numParam = varGroup->vars.size();
181         wanted = 0;
182         mac = parent? parent->mac : 0;
183         fit = parent? parent->fit : NA_REAL;
184         est = new double[numParam];
185         infoDefinite = NA_LOGICAL;
186         infoCondNum = NA_REAL;
187         infoA = NULL;
188         infoB = NULL;
189         stderrs = NULL;
190         inform = INFORM_UNINITIALIZED;
191         iterations = 0;
192
193         hess.resize(numParam, numParam);
194         ihess.resize(numParam, numParam);
195         clearHessian();
196 }
197
198 void FitContext::clearHessian()
199 {
200         for (size_t bx=0; bx < allBlocks.size(); ++bx) {
201                 delete allBlocks[bx];
202         }
203
204         allBlocks.clear();
205         haveSparseHess = false;
206         haveSparseIHess = false;
207         haveDenseHess = false;
208         haveDenseIHess = false;
209 }
210
211 void FitContext::allocStderrs()
212 {
213         if (stderrs) return;
214
215         size_t numParam = varGroup->vars.size();
216         stderrs = new double[numParam];
217
218         for (size_t px=0; px < numParam; ++px) {
219                 stderrs[px] = NA_REAL;
220         }
221 }
222
223 FitContext::FitContext(std::vector<double> &startingValues)
224 {
225         parent = NULL;
226         varGroup = Global->freeGroup[FREEVARGROUP_ALL];
227         init();
228
229         size_t numParam = varGroup->vars.size();
230         if (startingValues.size() != numParam) {
231                 Rf_error("Got %d starting values for %d parameters",
232                       startingValues.size(), numParam);
233         }
234         memcpy(est, startingValues.data(), sizeof(double) * numParam);
235 }
236
237 FitContext::FitContext(FitContext *parent, FreeVarGroup *varGroup)
238 {
239         this->parent = parent;
240         this->varGroup = varGroup;
241         init();
242
243         FreeVarGroup *src = parent->varGroup;
244         FreeVarGroup *dest = varGroup;
245         size_t dvars = varGroup->vars.size();
246         if (dvars == 0) return;
247         mapToParent.resize(dvars);
248
249         size_t d1 = 0;
250         for (size_t s1=0; s1 < src->vars.size(); ++s1) {
251                 if (src->vars[s1] != dest->vars[d1]) continue;
252                 mapToParent[d1] = s1;
253                 est[d1] = parent->est[s1];
254                 if (++d1 == dvars) break;
255         }
256         if (d1 != dvars) Rf_error("Parent free parameter group (id=%d) is not a superset of %d",
257                                src->id[0], dest->id[0]);
258
259         wanted = parent->wanted;
260         infoDefinite = parent->infoDefinite;
261         infoCondNum = parent->infoCondNum;
262         IterationError = parent->IterationError;
263         iterations = parent->iterations;
264 }
265
266 void FitContext::copyParamToModel(omxMatrix *mat)
267 { copyParamToModel(mat->currentState); }
268
269 void FitContext::copyParamToModel(omxMatrix *mat, double *at)
270 { copyParamToModel(mat->currentState, at); }
271
272 void FitContext::updateParent()
273 {
274         FreeVarGroup *src = varGroup;
275         FreeVarGroup *dest = parent->varGroup;
276         size_t svars = varGroup->vars.size();
277
278         parent->wanted |= wanted;
279         parent->fit = fit;
280         parent->IterationError = IterationError;
281         parent->mac = mac;
282         parent->infoDefinite = infoDefinite;
283         parent->infoCondNum = infoCondNum;
284         parent->iterations = iterations;
285
286         // rewrite using mapToParent TODO
287
288         if (svars > 0) {
289                 size_t s1 = 0;
290                 for (size_t d1=0; d1 < dest->vars.size(); ++d1) {
291                         if (dest->vars[d1] != src->vars[s1]) continue;
292                         parent->est[d1] = est[s1];
293                         if (++s1 == svars) break;
294                 }
295                 if (stderrs) {
296                         parent->allocStderrs();
297                         for (size_t s1=0; s1 < src->vars.size(); ++s1) {
298                                 parent->stderrs[mapToParent[s1]] = stderrs[s1];
299                         }
300                 }
301         }
302         
303         // pda(est, 1, svars);
304         // pda(parent->est, 1, dvars);
305 }
306
307 void FitContext::updateParentAndFree()
308 {
309         updateParent();
310         delete this;
311 }
312
313 void FitContext::log(int what)
314 {
315         size_t count = varGroup->vars.size();
316         std::string buf;
317         if (what & FF_COMPUTE_MAXABSCHANGE) buf += string_snprintf("MAC: %.5f\n", mac);
318         if (what & FF_COMPUTE_FIT) buf += string_snprintf("fit: %.5f (scale %f)\n", fit, Global->llScale);
319         if (what & FF_COMPUTE_ESTIMATE) {
320                 buf += string_snprintf("est %d: c(", (int) count);
321                 for (size_t vx=0; vx < count; ++vx) {
322                         buf += string_snprintf("%.5f", est[vx]);
323                         if (vx < count - 1) buf += ", ";
324                 }
325                 buf += ")\n";
326         }
327         mxLogBig(buf);
328 }
329
330 void FitContext::recordIterationError(const char* msg, ...)
331 {
332         // could record the parameter vector here for a better error message TODO
333         va_list ap;
334         va_start(ap, msg);
335         std::string str = string_vsnprintf(msg, ap);
336         va_end(ap);
337
338 #pragma omp critical(IterationError)
339         IterationError = str;
340 }
341
342 std::string FitContext::getIterationError()
343 {
344         return IterationError;
345 }
346
347 static void _fixSymmetry(const char *name, double *mat, size_t numParam, bool force)
348 {
349         for (size_t h1=1; h1 < numParam; h1++) {
350                 for (size_t h2=0; h2 < h1; h2++) {
351                         if (!force && mat[h2 * numParam + h1] != 0) {
352                                 omxRaiseErrorf("%s is not upper triangular", name);
353                                 break;
354                         }
355                         mat[h2 * numParam + h1] = mat[h1 * numParam + h2];
356                 }
357         }
358 }
359
360 static void omxRepopulateRFitFunction(omxFitFunction* oo, double* x, int n)
361 {
362         omxRFitFunction* rFitFunction = (omxRFitFunction*)oo->argStruct;
363
364         SEXP theCall, estimate;
365
366         Rf_protect(estimate = Rf_allocVector(REALSXP, n));
367         double *est = REAL(estimate);
368         for(int i = 0; i < n ; i++) {
369                 est[i] = x[i];
370         }
371
372         Rf_protect(theCall = Rf_allocVector(LANGSXP, 4));
373
374         SETCAR(theCall, Rf_install("imxUpdateModelValues"));
375         SETCADR(theCall, rFitFunction->model);
376         SETCADDR(theCall, rFitFunction->flatModel);
377         SETCADDDR(theCall, estimate);
378
379         R_Reprotect(rFitFunction->model = Rf_eval(theCall, R_GlobalEnv), rFitFunction->modelIndex);
380
381         Rf_unprotect(2); // theCall, estimate
382 }
383
384 void FitContext::copyParamToModel(omxState* os)
385 {
386         copyParamToModel(os, est);
387 }
388
389 void FitContext::copyParamToModel(omxState* os, double *at)
390 {
391         size_t numParam = varGroup->vars.size();
392
393         if(numParam == 0) return;
394
395         // Confidence Intervals & Hessian Calculation probe the parameter space
396         // near the best estimate. If stale, we need to restore the best estimate
397         // before returning results to the user.
398         os->stale = at != est;
399
400         os->computeCount++;
401
402         if(OMX_VERBOSE) {
403                 std::string buf;
404                 buf += string_snprintf("Call: %d (%d) ", iterations, (int) os->computeCount);
405                 buf += ("Estimates: [");
406                 for(size_t k = 0; k < numParam; k++) {
407                         buf += string_snprintf(" %f", at[k]);
408                 }
409                 buf += ("]\n");
410                 mxLogBig(buf);
411         }
412
413         for(size_t k = 0; k < numParam; k++) {
414                 omxFreeVar* freeVar = varGroup->vars[k];
415                 for(size_t l = 0; l < freeVar->locations.size(); l++) {
416                         omxFreeVarLocation *loc = &freeVar->locations[l];
417                         omxMatrix *matrix = os->matrixList[loc->matrix];
418                         int row = loc->row;
419                         int col = loc->col;
420                         omxSetMatrixElement(matrix, row, col, at[k]);
421                         if(OMX_DEBUG) {
422                                 mxLog("Setting location (%d, %d) of matrix %d to value %f for var %d",
423                                         row, col, loc->matrix, at[k], (int) k);
424                         }
425                 }
426         }
427
428         if (RFitFunction) omxRepopulateRFitFunction(RFitFunction, at, numParam);
429
430         varGroup->markDirty(os);
431
432         if (os->childList.size() == 0) return;
433
434         for(size_t i = 0; i < os->childList.size(); i++) {
435                 copyParamToModel(os->childList[i], at);
436         }
437 }
438
439 double *FitContext::take(int want)
440 {
441         if (!(want & (wanted | FF_COMPUTE_ESTIMATE))) {
442                 Rf_error("Attempt to take %d but not available", want);
443         }
444
445         double *ret = NULL;
446         switch(want) {
447         case FF_COMPUTE_ESTIMATE:
448                 ret = est;
449                 est = NULL;
450                 break;
451         default:
452                 Rf_error("Taking of %d is not implemented", want);
453         }
454         if (!ret) Rf_error("Attempt to take %d, already taken", want);
455         return ret;
456 }
457
458 // Rethink this whole design TODO
459 // If we compute things blockwise then most of this can go away?
460 void FitContext::preInfo()
461 {
462         size_t npsq = numParam * numParam;
463
464         if (!infoA) infoA = new double[npsq];
465         if (!infoB) infoB = new double[npsq];
466
467         switch (infoMethod) {
468         case INFO_METHOD_SANDWICH:
469         case INFO_METHOD_MEAT:
470                 OMXZERO(infoB, npsq);
471         case INFO_METHOD_BREAD:
472                 OMXZERO(infoA, npsq);
473                 break;
474         case INFO_METHOD_HESSIAN:
475                 clearHessian();
476                 break;
477         default:
478                 Rf_error("Unknown information matrix estimation method %d", infoMethod);
479         }
480 }
481
482 void FitContext::postInfo()
483 {
484         size_t numParam = varGroup->vars.size();
485         switch (infoMethod) {
486         case INFO_METHOD_SANDWICH:{
487                 // move into FCDeriv TODO
488                 omxBuffer<double> work(numParam * numParam);
489                 Matrix amat(infoA, numParam, numParam);
490                 InvertSymmetricIndef(amat, 'U');
491                 _fixSymmetry("InfoB", infoB, numParam, false);
492                 Matrix bmat(infoB, numParam, numParam);
493                 Matrix wmat(work.data(), numParam, numParam);
494                 Matrix hmat(getDenseIHessUninitialized(), numParam, numParam);
495                 SymMatrixMultiply('L', 'U', 1, 0, amat, bmat, wmat);
496                 SymMatrixMultiply('R', 'U', 1, 0, amat, wmat, hmat);
497                 wanted |= FF_COMPUTE_IHESSIAN;
498                 break;}
499         case INFO_METHOD_MEAT:{
500                 memcpy(getDenseHessUninitialized(), infoB, sizeof(double) * numParam * numParam); // avoid copy TODO
501                 wanted |= FF_COMPUTE_HESSIAN;
502                 break;}
503         case INFO_METHOD_BREAD:{
504                 memcpy(getDenseHessUninitialized(), infoA, sizeof(double) * numParam * numParam); // avoid copy TODO
505                 wanted |= FF_COMPUTE_HESSIAN;
506                 break;}
507         case INFO_METHOD_HESSIAN:
508                 if (Global->llScale > 0) negateHessian();
509                 wanted |= FF_COMPUTE_HESSIAN;
510                 break;
511         default:
512                 Rf_error("Unknown information matrix estimation method %d", infoMethod);
513         }
514 }
515
516 FitContext::~FitContext()
517 {
518         clearHessian();
519         if (est) delete [] est;
520         if (stderrs) delete [] stderrs;
521         if (infoA) delete [] infoA;
522         if (infoB) delete [] infoB;
523 }
524
525 omxFitFunction *FitContext::RFitFunction = NULL;
526
527 void FitContext::setRFitFunction(omxFitFunction *rff)
528 {
529         if (rff) {
530                 Global->numThreads = 1;
531                 if (RFitFunction) {
532                         Rf_error("You can only create 1 MxRFitFunction per independent model");
533                 }
534         }
535         RFitFunction = rff;
536 }
537
538 Ramsay1975::Ramsay1975(FitContext *fc, const char *flavor, int verbose, double minCaution) :
539         fc(fc), flavor(flavor), verbose(verbose), minCaution(minCaution)
540 {
541         if (!flavor) Rf_error("Ramsay: flavor cannot be NULL");
542
543         maxCaution = 0.0;
544         boundsHit = 0;
545         caution = 0;
546         highWatermark = std::max(0.5, caution);  // arbitrary guess
547
548         numParam = fc->varGroup->vars.size();
549         for (size_t px=0; px < numParam; ++px) {
550                 if (strcmp(fc->flavor[px], flavor) != 0) continue;
551                 vars.push_back(px);
552         }
553
554         prevAdj1.assign(numParam, 0);
555         prevAdj2.resize(numParam);
556         prevEst.resize(numParam);
557         memcpy(prevEst.data(), fc->est, sizeof(double) * numParam);
558
559         if (verbose >= 2) {
560                 mxLog("Ramsay[%10s]: %d parameters, caution %f, min caution %f",
561                       flavor, (int)vars.size(), caution, minCaution);
562         }
563 }
564
565 void Ramsay1975::recordEstimate(int px, double newEst)
566 {
567         omxFreeVar *fv = fc->varGroup->vars[px];
568         bool hitBound=false;
569         double param = newEst;
570         if (param < fv->lbound) {
571                 hitBound=true;
572                 param = prevEst[px] - (prevEst[px] - fv->lbound) / 2;
573         }
574         if (param > fv->ubound) {
575                 hitBound=true;
576                 param = prevEst[px] + (fv->ubound - prevEst[px]) / 2;
577         }
578         boundsHit += hitBound;
579         
580         prevAdj2[px] = prevAdj1[px];
581         prevAdj1[px] = param - prevEst[px];
582         
583         if (verbose >= 4) {
584                 std::string buf;
585                 buf += string_snprintf("Ramsay[%10s]: %d~%s %.4f -> %.4f", flavor, px, fv->name, prevEst[px], param);
586                 if (hitBound) {
587                         buf += string_snprintf(" wanted %.4f but hit bound", newEst);
588                 }
589                 if (prevAdj1[px] * prevAdj2[px] < 0) {
590                         buf += " *OSC*";
591                 }
592                 buf += "\n";
593                 mxLogBig(buf);
594         }
595
596         fc->est[px] = param;
597         prevEst[px] = param;
598 }
599
600 void Ramsay1975::apply()
601 {
602         for (size_t px=0; px < vars.size(); ++px) {
603                 int vx = vars[px];
604                 recordEstimate(vx, (1 - caution) * fc->est[vx] + caution * prevEst[vx]);
605         }
606 }
607
608 void Ramsay1975::recalibrate(bool *restart)
609 {
610         if (vars.size() == 0) return;
611
612         double normPrevAdj2 = 0;
613         double normAdjDiff = 0;
614         std::vector<double> adjDiff(numParam);
615
616         // The choice of norm is also arbitrary. Other norms might work better.
617         for (size_t vx=0; vx < vars.size(); ++vx) {
618                 int px = vars[vx];
619                 adjDiff[px] = prevAdj1[px] - prevAdj2[px];
620                 normPrevAdj2 += prevAdj2[px] * prevAdj2[px];
621         }
622
623         for (size_t vx=0; vx < vars.size(); ++vx) {
624                 int px = vars[vx];
625                 normAdjDiff += adjDiff[px] * adjDiff[px];
626         }
627         if (normAdjDiff == 0) {
628                 return;
629                 //Rf_error("Ramsay: no free variables of flavor %d", flavor);
630         }
631
632         double ratio = sqrt(normPrevAdj2 / normAdjDiff);
633         //if (verbose >= 3) mxLog("Ramsay[%d]: sqrt(%.5f/%.5f) = %.5f",
634         // flavor, normPrevAdj2, normAdjDiff, ratio);
635
636         double newCaution = 1 - (1-caution) * ratio / (1+boundsHit);
637         if (newCaution > .95) newCaution = .95;  // arbitrary guess
638         if (newCaution < 0) newCaution /= 2;     // don't get overconfident
639         if (newCaution < minCaution) newCaution = minCaution;
640         if (newCaution < caution) {
641                 caution = newCaution/3 + 2*caution/3;  // don't speed up too fast, arbitrary ratio
642         } else {
643                 caution = newCaution;
644         }
645         maxCaution = std::max(maxCaution, caution);
646         goingWild = false;
647         if (caution < highWatermark || (normPrevAdj2 < 1e-3 && normAdjDiff < 1e-3)) {
648                 if (verbose >= 3) mxLog("Ramsay[%10s]: %.2f caution", flavor, caution);
649         } else {
650                 if (verbose >= 3) {
651                         mxLog("Ramsay[%10s]: caution %.2f > %.2f, extreme oscillation, restart recommended",
652                               flavor, caution, highWatermark);
653                 }
654                 *restart = TRUE;
655                 goingWild = true;
656         }
657         highWatermark += .02; // arbitrary guess
658         boundsHit = 0;
659 }
660
661 void Ramsay1975::restart(bool myFault)
662 {
663         memcpy(prevEst.data(), fc->est, sizeof(double) * numParam);
664         prevAdj1.assign(numParam, 0);
665         prevAdj2.assign(numParam, 0);
666         myFault |= goingWild;
667         if (myFault) {
668                 highWatermark = 1 - (1 - highWatermark) * .5; // arbitrary guess
669                 caution = std::max(caution, highWatermark);   // arbitrary guess
670                 maxCaution = std::max(maxCaution, caution);
671                 highWatermark = caution;
672         }
673         if (vars.size() && verbose >= 3) {
674                 mxLog("Ramsay[%10s]: restart%s with %.2f caution %.2f highWatermark",
675                       flavor, myFault? " (my fault)":"", caution, highWatermark);
676         }
677 }
678
679 omxCompute::omxCompute()
680 {
681         varGroup = NULL;
682 }
683
684 void omxCompute::collectResultsHelper(FitContext *fc, std::vector< omxCompute* > &clist,
685                                       LocalComputeResult *lcr, MxRList *out)
686 {
687         for (std::vector< omxCompute* >::iterator it = clist.begin(); it != clist.end(); ++it) {
688                 omxCompute *c1 = *it;
689                 c1->collectResults(fc, lcr, out);
690         }
691 }
692
693 void omxCompute::collectResults(FitContext *fc, LocalComputeResult *lcr, MxRList *out)
694 {
695         MxRList *slots = new MxRList();
696         reportResults(fc, slots, out);
697         if (slots->size()) {
698                 lcr->push_back(std::make_pair(computeId, slots));
699         } else {
700                 delete slots;
701         }
702 }
703
704 omxCompute::~omxCompute()
705 {}
706
707 void omxCompute::initFromFrontend(SEXP rObj)
708 {
709         SEXP slotValue;
710         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("id")));
711         if (Rf_length(slotValue) != 1) Rf_error("MxCompute has no ID");
712
713         computeId = INTEGER(slotValue)[0];
714         varGroup = Global->findVarGroup(computeId);
715
716         if (!varGroup) {
717                 Rf_protect(slotValue = R_do_slot(rObj, Rf_install("freeSet")));
718                 if (Rf_length(slotValue) == 0) {
719                         varGroup = Global->findVarGroup(FREEVARGROUP_NONE);
720                 } else if (strcmp(CHAR(STRING_ELT(slotValue, 0)), ".")==0) {
721                         varGroup = Global->freeGroup[FREEVARGROUP_ALL];
722                 } else {
723                         Rf_warning("MxCompute ID %d references matrix '%s' in its freeSet "
724                                 "but this matrix contains no free parameters",
725                                 computeId, CHAR(STRING_ELT(slotValue, 0)));
726                         varGroup = Global->findVarGroup(FREEVARGROUP_NONE);
727                 }
728         }
729         if (OMX_DEBUG) {
730                 mxLog("MxCompute id %d assigned to var group %d", computeId, varGroup->id[0]);
731         }
732 }
733
734 void omxCompute::compute(FitContext *fc)
735 {
736         ComputeInform origInform = fc->inform;
737         FitContext *narrow = fc;
738         if (fc->varGroup != varGroup) narrow = new FitContext(fc, varGroup);
739         narrow->inform = INFORM_UNINITIALIZED;
740         computeImpl(narrow);
741         fc->inform = std::max(origInform, narrow->inform);
742         if (fc->varGroup != varGroup) narrow->updateParentAndFree();
743         Global->checkpointMessage(fc, fc->est, "%s", name);
744 }
745
746 class ComputeContainer : public omxCompute {
747         typedef omxCompute super;
748 protected:
749         std::vector< omxCompute* > clist;
750 public:
751         virtual void collectResults(FitContext *fc, LocalComputeResult *lcr, MxRList *out);
752 };
753
754 void ComputeContainer::collectResults(FitContext *fc, LocalComputeResult *lcr, MxRList *out)
755 {
756         super::collectResults(fc, lcr, out);
757         collectResultsHelper(fc, clist, lcr, out);
758 }
759
760 class omxComputeSequence : public ComputeContainer {
761         typedef ComputeContainer super;
762
763  public:
764         virtual void initFromFrontend(SEXP rObj);
765         virtual void computeImpl(FitContext *fc);
766         virtual ~omxComputeSequence();
767 };
768
769 class omxComputeIterate : public ComputeContainer {
770         typedef ComputeContainer super;
771         int maxIter;
772         double tolerance;
773         int verbose;
774
775  public:
776         virtual void initFromFrontend(SEXP rObj);
777         virtual void computeImpl(FitContext *fc);
778         virtual ~omxComputeIterate();
779 };
780
781 class omxComputeOnce : public omxCompute {
782         typedef omxCompute super;
783         std::vector< omxMatrix* > algebras;
784         std::vector< omxExpectation* > expectations;
785         std::vector< const char* > predict;
786         const char *how;
787         int verbose;
788         bool mac;
789         bool starting;
790         bool fit;
791         bool gradient;
792         bool hessian;
793         bool ihessian;
794         bool infoMat;
795         enum ComputeInfoMethod infoMethod;
796         bool hgprod;
797         bool isBestFit; // for backward compatibility
798
799  public:
800         virtual void initFromFrontend(SEXP rObj);
801         virtual void computeImpl(FitContext *fc);
802         virtual void reportResults(FitContext *fc, MxRList *slots, MxRList *out);
803 };
804
805 class ComputeEM : public omxCompute {
806         typedef omxCompute super;
807         std::vector< omxExpectation* > expectations;
808         const char *predict;
809         omxCompute *fit1;  // maybe rename to mstep TODO
810         omxMatrix *fit3;   // rename to observedFit
811         int EMcycles;
812         int maxIter;
813         int mstepIter;
814         int totalMstepIter;
815         double tolerance;
816         double semTolerance;
817         int verbose;
818         bool useRamsay;
819         bool information;
820         std::vector< omxMatrix * > semFitFunction;
821         enum ComputeInfoMethod infoMethod;
822         enum SEMMethod { ClassicSEM, TianSEM, GridSEM, AgileSEM } semMethod;
823         double *semMethodData;
824         int semMethodLen;
825         bool semDebug;
826         bool semFixSymmetry;
827         bool semForcePD;
828         int agileMaxIter;
829         SEXP rateMatrix; //debug
830         SEXP inputInfoMatrix; //debug
831         SEXP outputInfoMatrix; //debug
832         SEXP origEigenvalues; //debug
833         std::vector<Ramsay1975*> ramsay;
834         double noiseTarget;
835         double noiseTolerance;
836         std::vector<double*> estHistory;
837         std::vector<double> probeOffset;
838         std::vector<double> diffWork;
839         std::vector<int> paramHistLen;
840         std::vector<double> optimum;
841         double bestFit;
842         static const double MIDDLE_START;
843         static const double MIDDLE_END;
844         size_t maxHistLen;
845         int semProbeCount;
846
847         void setExpectationPrediction(const char *context);
848         void probeEM(FitContext *fc, int vx, double offset, std::vector<double> *rijWork);
849         void recordDiff(FitContext *fc, int v1, std::vector<double> &rijWork,
850                         double *stdDiff, bool *mengOK);
851
852  public:
853         virtual void initFromFrontend(SEXP rObj);
854         virtual void computeImpl(FitContext *fc);
855         virtual void collectResults(FitContext *fc, LocalComputeResult *lcr, MxRList *out);
856         virtual void reportResults(FitContext *fc, MxRList *slots, MxRList *out);
857         virtual ~ComputeEM();
858 };
859
860 const double ComputeEM::MIDDLE_START = 0.105360515657826281366; // -log(.9) constexpr
861 const double ComputeEM::MIDDLE_END = 0.001000500333583534363566; // -log(.999) constexpr
862
863 class ComputeStandardError : public omxCompute {
864         typedef omxCompute super;
865  public:
866         virtual void reportResults(FitContext *fc, MxRList *slots, MxRList *out);
867 };
868
869 class ComputeHessianQuality : public omxCompute {
870         typedef omxCompute super;
871  public:
872         virtual void reportResults(FitContext *fc, MxRList *slots, MxRList *out);
873 };
874
875 class ComputeReportDeriv : public omxCompute {
876         typedef omxCompute super;
877  public:
878         virtual void reportResults(FitContext *fc, MxRList *slots, MxRList *out);
879 };
880
881 static class omxCompute *newComputeSequence()
882 { return new omxComputeSequence(); }
883
884 static class omxCompute *newComputeIterate()
885 { return new omxComputeIterate(); }
886
887 static class omxCompute *newComputeOnce()
888 { return new omxComputeOnce(); }
889
890 static class omxCompute *newComputeEM()
891 { return new ComputeEM(); }
892
893 static class omxCompute *newComputeStandardError()
894 { return new ComputeStandardError(); }
895
896 static class omxCompute *newComputeHessianQuality()
897 { return new ComputeHessianQuality(); }
898
899 static class omxCompute *newComputeReportDeriv()
900 { return new ComputeReportDeriv(); }
901
902 struct omxComputeTableEntry {
903         char name[32];
904         omxCompute *(*ctor)();
905 };
906
907 static const struct omxComputeTableEntry omxComputeTable[] = {
908         {"MxComputeNumericDeriv", &newComputeNumericDeriv},
909         {"MxComputeGradientDescent", &newComputeGradientDescent},
910         {"MxComputeSequence", &newComputeSequence },
911         {"MxComputeIterate", &newComputeIterate },
912         {"MxComputeOnce", &newComputeOnce },
913         {"MxComputeNewtonRaphson", &newComputeNewtonRaphson},
914         {"MxComputeEM", &newComputeEM },
915         {"MxComputeStandardError", &newComputeStandardError},
916         {"MxComputeHessianQuality", &newComputeHessianQuality},
917         {"MxComputeReportDeriv", &newComputeReportDeriv},
918         {"MxComputeConfidenceInterval", &newComputeConfidenceInterval}
919 };
920
921 omxCompute *omxNewCompute(omxState* os, const char *type)
922 {
923         omxCompute *got = NULL;
924
925         for (size_t fx=0; fx < OMX_STATIC_ARRAY_SIZE(omxComputeTable); fx++) {
926                 const struct omxComputeTableEntry *entry = omxComputeTable + fx;
927                 if(strcmp(type, entry->name) == 0) {
928                         got = entry->ctor();
929                         got->name = entry->name;
930                         break;
931                 }
932         }
933
934         if (!got) Rf_error("Compute %s is not implemented", type);
935
936         return got;
937 }
938
939 void omxComputeSequence::initFromFrontend(SEXP rObj)
940 {
941         super::initFromFrontend(rObj);
942
943         SEXP slotValue;
944         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("steps")));
945
946         for (int cx = 0; cx < Rf_length(slotValue); cx++) {
947                 SEXP step = VECTOR_ELT(slotValue, cx);
948                 SEXP s4class;
949                 Rf_protect(s4class = STRING_ELT(Rf_getAttrib(step, Rf_install("class")), 0));
950                 omxCompute *compute = omxNewCompute(globalState, CHAR(s4class));
951                 compute->initFromFrontend(step);
952                 if (isErrorRaised(globalState)) break;
953                 clist.push_back(compute);
954         }
955 }
956
957 void omxComputeSequence::computeImpl(FitContext *fc)
958 {
959         for (size_t cx=0; cx < clist.size(); ++cx) {
960                 clist[cx]->compute(fc);
961                 if (isErrorRaised(globalState)) break;
962         }
963 }
964
965 omxComputeSequence::~omxComputeSequence()
966 {
967         for (size_t cx=0; cx < clist.size(); ++cx) {
968                 delete clist[cx];
969         }
970 }
971
972 void omxComputeIterate::initFromFrontend(SEXP rObj)
973 {
974         SEXP slotValue;
975
976         super::initFromFrontend(rObj);
977
978         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("maxIter")));
979         maxIter = INTEGER(slotValue)[0];
980
981         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("tolerance")));
982         tolerance = REAL(slotValue)[0];
983         if (tolerance <= 0) Rf_error("tolerance must be positive");
984
985         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("steps")));
986
987         for (int cx = 0; cx < Rf_length(slotValue); cx++) {
988                 SEXP step = VECTOR_ELT(slotValue, cx);
989                 SEXP s4class;
990                 Rf_protect(s4class = STRING_ELT(Rf_getAttrib(step, Rf_install("class")), 0));
991                 omxCompute *compute = omxNewCompute(globalState, CHAR(s4class));
992                 compute->initFromFrontend(step);
993                 if (isErrorRaised(globalState)) break;
994                 clist.push_back(compute);
995         }
996
997         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("verbose")));
998         verbose = Rf_asInteger(slotValue);
999 }
1000
1001 void omxComputeIterate::computeImpl(FitContext *fc)
1002 {
1003         int iter = 0;
1004         double prevFit = 0;
1005         double mac = tolerance * 10;
1006         while (1) {
1007                 ++fc->iterations;
1008                 for (size_t cx=0; cx < clist.size(); ++cx) {
1009                         clist[cx]->compute(fc);
1010                         if (isErrorRaised(globalState)) break;
1011                 }
1012                 if (fc->wanted & FF_COMPUTE_MAXABSCHANGE) {
1013                         if (fc->mac < 0) {
1014                                 Rf_warning("MAC estimated at %.4f; something is wrong", fc->mac);
1015                                 break;
1016                         } else {
1017                                 mac = fc->mac;
1018                                 if (verbose) mxLog("ComputeIterate: mac %.9g", mac);
1019                         }
1020                 }
1021                 if (fc->wanted & FF_COMPUTE_FIT) {
1022                         if (fc->fit == 0) {
1023                                 Rf_warning("Fit estimated at 0; something is wrong");
1024                                 break;
1025                         }
1026                         if (prevFit != 0) {
1027                                 double change = prevFit - fc->fit;
1028                                 if (verbose) mxLog("ComputeIterate: fit %.9g change %.9g", fc->fit, change);
1029                                 mac = fabs(change);
1030                         } else {
1031                                 if (verbose) mxLog("ComputeIterate: initial fit %.9g", fc->fit);
1032                         }
1033                         prevFit = fc->fit;
1034                 }
1035                 if (!(fc->wanted & (FF_COMPUTE_MAXABSCHANGE | FF_COMPUTE_FIT))) {
1036                         omxRaiseErrorf("ComputeIterate: neither MAC nor fit available");
1037                 }
1038                 if (isErrorRaised(globalState) || ++iter > maxIter || mac < tolerance) break;
1039         }
1040 }
1041
1042 omxComputeIterate::~omxComputeIterate()
1043 {
1044         for (size_t cx=0; cx < clist.size(); ++cx) {
1045                 delete clist[cx];
1046         }
1047 }
1048
1049 void ComputeEM::initFromFrontend(SEXP rObj)
1050 {
1051         SEXP slotValue;
1052         SEXP s4class;
1053
1054         super::initFromFrontend(rObj);
1055
1056         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("expectation")));
1057         for (int wx=0; wx < Rf_length(slotValue); ++wx) {
1058                 int objNum = INTEGER(slotValue)[wx];
1059                 omxExpectation *expectation = globalState->expectationList[objNum];
1060                 setFreeVarGroup(expectation, varGroup);
1061                 omxCompleteExpectation(expectation);
1062                 expectations.push_back(expectation);
1063         }
1064
1065         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("predict")));
1066         {
1067                 // Should accept a vector here TODO
1068                 if (Rf_length(slotValue) != 1) Rf_error("Not implemented");
1069                 SEXP elem;
1070                 Rf_protect(elem = STRING_ELT(slotValue, 0));
1071                 predict = CHAR(elem);
1072         }
1073
1074         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("mstep")));
1075         Rf_protect(s4class = STRING_ELT(Rf_getAttrib(slotValue, Rf_install("class")), 0));
1076         fit1 = omxNewCompute(globalState, CHAR(s4class));
1077         fit1->initFromFrontend(slotValue);
1078
1079         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("observedFit")));
1080         fit3 = globalState->algebraList[ INTEGER(slotValue)[0] ];
1081         if (fit3->fitFunction) {
1082                 setFreeVarGroup(fit3->fitFunction, varGroup);
1083                 omxCompleteFitFunction(fit3);
1084         }
1085
1086         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("maxIter")));
1087         maxIter = INTEGER(slotValue)[0];
1088
1089         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("tolerance")));
1090         tolerance = REAL(slotValue)[0];
1091         if (tolerance <= 0) Rf_error("tolerance must be positive");
1092
1093         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("verbose")));
1094         verbose = Rf_asInteger(slotValue);
1095
1096         useRamsay = false;
1097         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("accel")));
1098         const char *accelName = CHAR(STRING_ELT(slotValue, 0));
1099         if (strEQ(accelName, "ramsay1975")) {
1100                 useRamsay = true;
1101         } else if (STRING_ELT(slotValue, 0) == NA_STRING) {
1102                 // OK
1103         } else {
1104                 Rf_warning("%s: unknown acceleration method %s ignored", name, accelName);
1105         }
1106
1107         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("information")));
1108         const char *infoName = CHAR(STRING_ELT(slotValue, 0));
1109         information = false;
1110         if (STRING_ELT(slotValue, 0) == NA_STRING) {
1111                 // ok
1112         } else if (strEQ(infoName, "mr1991")) {
1113                 information = true;
1114         } else {
1115                 Rf_warning("%s: unknown information method %s ignored", name, infoName);
1116         }
1117
1118         if (information) {
1119                 infoMethod = INFO_METHOD_HESSIAN;
1120                 semMethod = AgileSEM;
1121                 agileMaxIter = 1;
1122                 semDebug = false;
1123                 semFixSymmetry = true;
1124                 semForcePD = false;
1125                 noiseTarget = exp(-5.0); //constexpr
1126                 noiseTolerance = exp(3.3); //constexpr
1127                 semTolerance = sqrt(tolerance);  // override needed?
1128
1129                 SEXP infoArgs, argNames;
1130                 Rf_protect(infoArgs = R_do_slot(rObj, Rf_install("infoArgs")));
1131                 Rf_protect(argNames = Rf_getAttrib(infoArgs, R_NamesSymbol));
1132
1133                 for (int ax=0; ax < Rf_length(infoArgs); ++ax) {
1134                         const char *key = R_CHAR(STRING_ELT(argNames, ax));
1135                         slotValue = VECTOR_ELT(infoArgs, ax);
1136                         if (strEQ(key, "fitfunction")) {
1137                                 for (int fx=0; fx < Rf_length(slotValue); ++fx) {
1138                                         omxMatrix *ff = globalState->algebraList[INTEGER(slotValue)[fx]];
1139                                         if (!ff->fitFunction) Rf_error("infoArgs$fitfunction is %s, not a fitfunction", ff->name);
1140                                         semFitFunction.push_back(ff);
1141                                 }
1142                         } else if (strEQ(key, "inputInfo")) {
1143                                 infoMethod = stringToInfoMethod(CHAR(slotValue));
1144                         } else if (strEQ(key, "semMethod")) {
1145                                 semMethodLen = Rf_length(slotValue);
1146                                 if (semMethodLen == 0) {
1147                                         semMethod = AgileSEM;
1148                                         semMethodData = NULL;
1149                                 } else {
1150                                         semMethodData = REAL(slotValue);
1151                                         if (semMethodLen > 1) {
1152                                                 semMethod = GridSEM;
1153                                         } else if (semMethodData[0] == 1) {
1154                                                 semMethod = ClassicSEM;
1155                                         } else if (semMethodData[0] == 2) {
1156                                                 semMethod = TianSEM;
1157                                         } else if (semMethodData[0] == 3) {
1158                                                 semMethod = AgileSEM;
1159                                         } else {
1160                                                 Rf_error("Unknown SEM method %f", semMethodData[0]);
1161                                         }
1162                                 }
1163                         } else if (strEQ(key, "agileMaxIter")) {
1164                                 agileMaxIter = INTEGER(slotValue)[0];
1165                         } else if (strEQ(key, "semDebug")) {
1166                                 semDebug = Rf_asLogical(slotValue);
1167                         } else if (strEQ(key, "semFixSymmetry")) {
1168                                 semFixSymmetry = Rf_asLogical(slotValue);
1169                         } else if (strEQ(key, "semForcePD")) {
1170                                 semForcePD = Rf_asLogical(slotValue);
1171                         } else if (strEQ(key, "noiseTarget")) {
1172                                 noiseTarget = REAL(slotValue)[0];
1173                                 if (noiseTarget <= 0) Rf_error("noiseTarget must be positive");
1174                         } else if (strEQ(key, "noiseTolerance")) {
1175                                 noiseTolerance = REAL(slotValue)[0];
1176                                 if (noiseTolerance < 1) Rf_error("noiseTolerance must be >=1");
1177                         } else {
1178                                 mxLog("%s: unknown key %s", name, key);
1179                         }
1180                 }
1181                 if (!semFixSymmetry && semForcePD) {
1182                         Rf_warning("%s: semFixSymmetry must be enabled for semForcePD", name);
1183                         semForcePD = false;
1184                 }
1185         }
1186
1187         inputInfoMatrix = NULL;
1188         outputInfoMatrix = NULL;
1189         rateMatrix = NULL;
1190         origEigenvalues = NULL;
1191 }
1192
1193 void ComputeEM::setExpectationPrediction(const char *context)
1194 {
1195         for (size_t wx=0; wx < expectations.size(); ++wx) {
1196                 omxExpectation *expectation = expectations[wx];
1197                 if (verbose >= 4) mxLog("ComputeEM: expectation[%d] %s predict %s", (int) wx, expectation->name, context);
1198                 omxExpectationCompute(expectation, context);
1199         }
1200 }
1201
1202 void ComputeEM::probeEM(FitContext *fc, int vx, double offset, std::vector<double> *rijWork)
1203 {
1204         const size_t freeVars = fc->varGroup->vars.size();
1205         const int base = paramHistLen[vx] * freeVars;
1206         probeOffset[vx * maxHistLen + paramHistLen[vx]] = offset;
1207         paramHistLen[vx] += 1;
1208
1209         memcpy(fc->est, optimum.data(), sizeof(double) * freeVars);
1210         fc->est[vx] += offset;
1211         fc->copyParamToModel(globalState);
1212
1213         setExpectationPrediction(predict);
1214         int informSave = fc->inform;  // not sure if we want to hide inform here TODO
1215         fit1->compute(fc);
1216         fc->inform = informSave;
1217         setExpectationPrediction("nothing");
1218
1219         if (verbose >= 3) mxLog("ComputeEM: probe %d of param %d offset %.6f",
1220                                 paramHistLen[vx], vx, offset);
1221
1222         for (size_t v1=0; v1 < freeVars; ++v1) {
1223                 double got = (fc->est[v1] - optimum[v1]) / offset;
1224                 (*rijWork)[base + v1] = got;
1225         }
1226         //pda(rij.data() + base, 1, freeVars);
1227         ++semProbeCount;
1228 }
1229
1230 void ComputeEM::recordDiff(FitContext *fc, int v1, std::vector<double> &rijWork,
1231                            double *stdDiff, bool *mengOK)
1232 {
1233         const size_t freeVars = fc->varGroup->vars.size();
1234         int h1 = paramHistLen[v1]-2;
1235         int h2 = paramHistLen[v1]-1;
1236         double *rij1 = rijWork.data() + h1 * freeVars;
1237         double *rij2 = rijWork.data() + h2 * freeVars;
1238         double diff = 0;
1239         *mengOK = true;
1240         for (size_t v2=0; v2 < freeVars; ++v2) {
1241                 double diff1 = fabs(rij1[v2] - rij2[v2]);
1242                 if (diff1 >= semTolerance) *mengOK = false;
1243                 diff += diff1;
1244         }
1245         double p1 = probeOffset[v1 * maxHistLen + h1];
1246         double p2 = probeOffset[v1 * maxHistLen + h2];
1247         double dist = fabs(p1 - p2);
1248         if (dist < tolerance/4) Rf_error("SEM: invalid probe offset distance %.9f", dist);
1249         *stdDiff = diff / (freeVars * dist);
1250         diffWork[v1 * maxHistLen + h1] = *stdDiff;
1251         if (verbose >= 2) mxLog("ComputeEM: (%f,%f) mengOK %d diff %f stdDiff %f",
1252                                 p1, p2, *mengOK, diff / freeVars, *stdDiff);
1253 }
1254
1255 void ComputeEM::computeImpl(FitContext *fc)
1256 {
1257         const double Scale = fabs(Global->llScale);
1258         double prevFit = 0;
1259         double mac = tolerance * 10;
1260         bool converged = false;
1261         const size_t freeVars = fc->varGroup->vars.size();
1262         bool in_middle = false;
1263         maxHistLen = 0;
1264         EMcycles = 0;
1265         semProbeCount = 0;
1266
1267         if (verbose >= 1) mxLog("ComputeEM: Welcome, tolerance=%g ramsay=%d info=%d",
1268                                 tolerance, useRamsay, information);
1269
1270         // Some evidence suggests that better performance is obtained when the
1271         // item parameters and latent distribution parameters are split into
1272         // separate Ramsay1975 groups with different minimum caution limits,
1273         //
1274         // ramsay.push_back(new Ramsay1975(fc, 1+int(ramsay.size()), 0, verbose, -1.25)); // M-step param
1275         // ramsay.push_back(new Ramsay1975(fc, 1+int(ramsay.size()), 0, verbose, -1));    // extra param
1276         //
1277         // I had this hardcoded for a while, but in making the API more generic,
1278         // I'm not sure how to allow specification of the Ramsay1975 grouping.
1279         // One possibility is list(flavor1=c("ItemParam"), flavor2=c("mean","cov"))
1280         // but this doesn't allow finer grain than matrix-wise assignment. The
1281         // other question is whether more Ramsay1975 groups really help or not.
1282         // nightly/ifa-cai2009.R actually got faster with 1 Ramsay group.
1283
1284         const char *flavor = "EM";
1285         fc->flavor.assign(freeVars, flavor);
1286         ramsay.push_back(new Ramsay1975(fc, flavor, verbose, -1.25));
1287
1288         while (1) {
1289                 if (verbose >= 4) mxLog("ComputeEM[%d]: E-step", EMcycles);
1290                 setExpectationPrediction(predict);
1291
1292                 {
1293                         if (verbose >= 4) mxLog("ComputeEM[%d]: M-step", EMcycles);
1294                         int informSave = fc->inform;
1295                         FitContext *fc1 = new FitContext(fc, fit1->varGroup);
1296                         int startIter = fc1->iterations;
1297                         fit1->compute(fc1);
1298                         mstepIter = fc1->iterations - startIter;
1299                         fc1->updateParentAndFree();
1300                         fc->inform = informSave;
1301                 }
1302
1303                 {
1304                         if (useRamsay) {
1305                                 bool wantRestart;
1306                                 if (EMcycles > 3 && EMcycles % 3 == 0) {
1307                                         for (size_t rx=0; rx < ramsay.size(); ++rx) {
1308                                                 ramsay[rx]->recalibrate(&wantRestart);
1309                                         }
1310                                 }
1311                                 for (size_t rx=0; rx < ramsay.size(); ++rx) {
1312                                         ramsay[rx]->apply();
1313                                 }
1314                         }
1315                         fc->copyParamToModel(globalState);
1316                         if (verbose >= 4) mxLog("ComputeEM[%d]: observed fit", EMcycles);
1317                         setExpectationPrediction("nothing");
1318                         Global->checkpointPrefit(fc, fc->est, false);
1319                         omxForceCompute(fit3);
1320                         fc->fit = fit3->data[0];
1321                         Global->checkpointPostfit(fc);
1322                 }
1323
1324                 totalMstepIter += mstepIter;
1325
1326                 if (!(fc->wanted & FF_COMPUTE_FIT)) {
1327                         omxRaiseErrorf("ComputeEM: fit not available");
1328                         break;
1329                 }
1330                 if (fc->fit == 0) {
1331                         omxRaiseErrorf("Fit estimated at 0; something is wrong");
1332                         break;
1333                 }
1334                 double change = 0;
1335                 if (prevFit != 0) {
1336                         change = prevFit - fc->fit;
1337                         if (verbose >= 2) mxLog("ComputeEM[%d]: msteps %d fit %.9g change %.9g",
1338                                                 EMcycles, mstepIter, fc->fit, change);
1339                         mac = fabs(change);
1340                         if (mac < MIDDLE_START * Scale) in_middle = true;
1341                         if (mac < MIDDLE_END * Scale) in_middle = false;
1342                 } else {
1343                         if (verbose >= 2) mxLog("ComputeEM: msteps %d initial fit %.9g",
1344                                                 mstepIter, fc->fit);
1345                 }
1346
1347                 prevFit = fc->fit;
1348                 converged = mac < tolerance;
1349                 ++fc->iterations;
1350                 if (isErrorRaised(globalState) || ++EMcycles > maxIter || converged) break;
1351
1352                 if (semMethod == ClassicSEM || ((semMethod == TianSEM || semMethod == AgileSEM) && in_middle)) {
1353                         double *estCopy = new double[freeVars];
1354                         memcpy(estCopy, fc->est, sizeof(double) * freeVars);
1355                         estHistory.push_back(estCopy);
1356                 }
1357         }
1358
1359         int wanted = FF_COMPUTE_FIT | FF_COMPUTE_BESTFIT | FF_COMPUTE_ESTIMATE;
1360         fc->wanted = wanted;
1361         bestFit = fc->fit;
1362         if (verbose >= 1) mxLog("ComputeEM: cycles %d/%d total mstep %d fit %f",
1363                                 EMcycles, maxIter, totalMstepIter, bestFit);
1364
1365         if (!converged || !information) return;
1366
1367         if (verbose >= 1) mxLog("ComputeEM: tolerance=%f semMethod=%d, semTolerance=%f ideal noise=[%f,%f]",
1368                                 tolerance, semMethod, semTolerance,
1369                                 noiseTarget/noiseTolerance, noiseTarget*noiseTolerance);
1370
1371         optimum.resize(freeVars);
1372         memcpy(optimum.data(), fc->est, sizeof(double) * freeVars);
1373
1374         if (semMethod == AgileSEM) {
1375                 maxHistLen = 2 + agileMaxIter * 2;
1376         } else if (semMethod == ClassicSEM || semMethod == TianSEM) {
1377                 maxHistLen = estHistory.size();
1378         } else {
1379                 maxHistLen = semMethodLen;
1380         }
1381
1382         probeOffset.resize(maxHistLen * freeVars);
1383         diffWork.resize(maxHistLen * freeVars);
1384         paramHistLen.assign(freeVars, 0);
1385         omxBuffer<double> rij(freeVars * freeVars);
1386
1387         size_t semConverged=0;
1388         for (size_t v1=0; v1 < freeVars; ++v1) {
1389                 std::vector<double> rijWork(freeVars * maxHistLen);
1390                 int pick = 0;
1391                 bool paramConverged = false;
1392                 if (semMethod == AgileSEM) {
1393                         const double stepSize = tolerance;
1394
1395                         double offset1 = tolerance * 50;
1396                         double sign = 1;
1397                         if (estHistory.size()) {
1398                                 int hpick = estHistory.size() /2;
1399                                 double popt = optimum[v1];
1400                                 sign = (popt < estHistory[hpick][v1])? 1 : -1;
1401                                 offset1 = fabs(estHistory[hpick][v1] - popt);
1402                                 if (offset1 < 10 * tolerance) offset1 = 10 * tolerance;
1403                                 if (offset1 > 1000 * tolerance) offset1 = 1000 * tolerance;
1404                         }
1405
1406                         probeEM(fc, v1, sign * offset1, &rijWork);
1407                         double offset2 = offset1 + stepSize;
1408                         probeEM(fc, v1, sign * offset2, &rijWork);
1409                         double diff;
1410                         bool mengOK;
1411                         recordDiff(fc, v1, rijWork, &diff, &mengOK);
1412                         double midOffset = (offset1 + offset2) / 2;
1413
1414                         int iter = 0;
1415                         omxBuffer<double> coefHist(agileMaxIter);
1416                         while (++iter <= agileMaxIter &&
1417                                !(noiseTarget/noiseTolerance < diff && diff < noiseTarget*noiseTolerance)) {
1418                                 coefHist[iter-1] = diff * midOffset * midOffset;
1419                                 double coef = 0;
1420                                 for (int cx=0; cx < iter; ++cx) coef += coefHist[cx];
1421                                 coef /= iter;
1422                                 if (verbose >= 4) mxLog("ComputeEM: agile iter[%d] coef=%.6g", iter, coef);
1423                                 offset1 = sqrt(coef/noiseTarget);
1424                                 probeEM(fc, v1, sign * offset1, &rijWork);
1425                                 if (iter < agileMaxIter || semDebug) {
1426                                         offset2 = offset1 + stepSize;
1427                                         probeEM(fc, v1, sign * offset2, &rijWork);
1428                                         midOffset = (offset1 + offset2) / 2;
1429                                         recordDiff(fc, v1, rijWork, &diff, &mengOK);
1430                                 }
1431                                 pick += 2;
1432                         }
1433                         paramConverged = true;
1434                 } else if (semMethod == ClassicSEM || semMethod == TianSEM) {
1435                         if (!estHistory.size()) {
1436                                 if (verbose >= 1) mxLog("ComputeEM: no history available;"
1437                                                         " Classic or Tian SEM require convergence history");
1438                                 return;
1439                         }
1440                         for (size_t hx=0; hx < estHistory.size(); ++hx) {
1441                                 double popt = optimum[v1];
1442                                 double offset1 = estHistory[hx][v1] - popt;
1443                                 if (paramHistLen[v1] && fabs(probeOffset[v1 * maxHistLen + paramHistLen[v1]-1] -
1444                                                              offset1) < tolerance) continue;
1445                                 if (fabs(offset1) < tolerance) continue;
1446                                 probeEM(fc, v1, offset1, &rijWork);
1447                                 if (hx == 0) continue;
1448                                 pick = hx;
1449                                 double diff;
1450                                 bool mengOK;
1451                                 recordDiff(fc, v1, rijWork, &diff, &mengOK);
1452                                 if (mengOK) {
1453                                         paramConverged = true;
1454                                         break;
1455                                 }
1456                         }
1457                 } else {
1458                         for (int hx=0; hx < semMethodLen; ++hx) {
1459                                 probeEM(fc, v1, semMethodData[hx], &rijWork);
1460                                 if (hx == 0) continue;
1461                                 double diff;
1462                                 bool mengOK;
1463                                 recordDiff(fc, v1, rijWork, &diff, &mengOK);
1464                         }
1465                         paramConverged = true;
1466                 }
1467
1468                 if (paramConverged) {
1469                         ++semConverged;
1470                         memcpy(rij.data() + v1 * freeVars, rijWork.data() + pick*freeVars, sizeof(double) * freeVars);
1471                         if (verbose >= 2) mxLog("ComputeEM: param %d converged in %d probes",
1472                                                 (int) v1, paramHistLen[v1]);
1473                 } else {
1474                         if (verbose >= 2) mxLog("ComputeEM: param %d failed to converge after %d probes",
1475                                                 (int) v1, paramHistLen[v1]);
1476                         break;
1477                 }
1478         }
1479
1480         if (verbose >= 1) {
1481                 if (semConverged == freeVars) {
1482                         mxLog("ComputeEM: %d probes used to estimate Hessian", semProbeCount);
1483                 } else {
1484                         mxLog("ComputeEM: %d probes used for SEM but failed to converge", semProbeCount);
1485                 }
1486         }
1487         if (semConverged < freeVars) return;
1488
1489         fc->fit = bestFit;
1490         memcpy(fc->est, optimum.data(), sizeof(double) * freeVars);
1491         fc->copyParamToModel(globalState);
1492
1493         if (semDebug) {
1494                 Rf_protect(rateMatrix = Rf_allocMatrix(REALSXP, freeVars, freeVars));
1495                 memcpy(REAL(rateMatrix), rij.data(), sizeof(double) * freeVars * freeVars);
1496         }
1497
1498         // rij = I-rij
1499         for (size_t v1=0; v1 < freeVars; ++v1) {
1500                 for (size_t v2=0; v2 < freeVars; ++v2) {
1501                         int cell = v1 * freeVars + v2;
1502                         double entry = rij[cell];
1503                         if (v1 == v2) entry = 1 - entry;
1504                         else entry = -entry;
1505                         rij[cell] = entry;
1506                 }
1507         }
1508
1509         //mxLog("rij symm");
1510         //pda(rij.data(), freeVars, freeVars);
1511
1512         setExpectationPrediction(predict);
1513         fc->wanted = 0;
1514         fc->infoMethod = infoMethod;
1515         fc->preInfo();
1516         for (size_t fx=0; fx < semFitFunction.size(); ++fx) {
1517                 omxFitFunctionCompute(semFitFunction[fx]->fitFunction, FF_COMPUTE_INFO, fc);
1518         }
1519         fc->postInfo();
1520         setExpectationPrediction("nothing");
1521
1522         Rf_protect(inputInfoMatrix = Rf_allocMatrix(REALSXP, freeVars, freeVars));
1523         double *hess = REAL(inputInfoMatrix);
1524         fc->copyDenseHess(hess);
1525
1526         Matrix rijMat(rij.data(), freeVars, freeVars);
1527         Matrix hessMat(hess, freeVars, freeVars);
1528         omxBuffer<double> infoBuf(freeVars * freeVars);
1529         Matrix infoMat(infoBuf.data(), freeVars, freeVars);
1530
1531         SymMatrixMultiply('L', 'U', 1, 0, hessMat, rijMat, infoMat);  // result not symmetric!
1532
1533         double *ihess = fc->getDenseIHessUninitialized();
1534         int singular;
1535         if (semFixSymmetry) {
1536                 MeanSymmetric(infoMat);
1537                 singular = InvertSymmetricIndef(infoMat, 'U');
1538                 memcpy(ihess, infoBuf.data(), sizeof(double) * freeVars * freeVars);
1539         } else {
1540                 Matrix ihessMat(ihess, freeVars, freeVars);
1541                 singular = MatrixSolve(infoMat, ihessMat, true);
1542         }
1543         if (singular) {
1544                 if (verbose >= 1) mxLog("ComputeEM: SEM Hessian is singular %d", singular);
1545                 return;
1546         }
1547
1548         if (semForcePD) {
1549                 double *oev = NULL;
1550                 if (semDebug) {
1551                         Rf_protect(origEigenvalues = Rf_allocVector(REALSXP, freeVars));
1552                         oev = REAL(origEigenvalues);
1553                 }
1554                 Matrix mat(ihess, freeVars, freeVars);
1555                 InplaceForcePosSemiDef(mat, oev, &fc->infoCondNum);
1556         }
1557
1558         if (semDebug) {
1559                 // ihess is always symmetric, this could be asymmetric
1560                 Rf_protect(outputInfoMatrix = Rf_allocMatrix(REALSXP, freeVars, freeVars));
1561                 memcpy(REAL(outputInfoMatrix), ihess, sizeof(double) * freeVars * freeVars);
1562         }
1563
1564         fc->wanted = wanted | FF_COMPUTE_IHESSIAN;
1565         //pda(fc->ihess, freeVars, freeVars);
1566 }
1567
1568 void ComputeEM::collectResults(FitContext *fc, LocalComputeResult *lcr, MxRList *out)
1569 {
1570         super::collectResults(fc, lcr, out);
1571
1572         std::vector< omxCompute* > clist(1);
1573         clist[0] = fit1;
1574
1575         collectResultsHelper(fc, clist, lcr, out);
1576 }
1577
1578 void ComputeEM::reportResults(FitContext *fc, MxRList *slots, MxRList *)
1579 {
1580         size_t numFree = fc->varGroup->vars.size();
1581         if (!numFree) return;
1582
1583         MxRList out;
1584         out.add("EMcycles", Rf_ScalarInteger(EMcycles));
1585         out.add("totalMstep", Rf_ScalarInteger(totalMstepIter));
1586         out.add("semProbeCount", Rf_ScalarInteger(semProbeCount));
1587         slots->add("output", out.asR());
1588
1589         if (semDebug) {
1590                 const int freeVars = (int) fc->varGroup->vars.size();
1591                 MxRList dbg;
1592
1593                 if (probeOffset.size()) {
1594                         SEXP Rpo;
1595                         Rf_protect(Rpo = Rf_allocMatrix(REALSXP, maxHistLen, freeVars));
1596                         memcpy(REAL(Rpo), probeOffset.data(), sizeof(double) * maxHistLen * freeVars);
1597                         dbg.add("probeOffset", Rpo);
1598                 }
1599
1600                 if (diffWork.size()) {
1601                         SEXP Rdiff;
1602                         Rf_protect(Rdiff = Rf_allocMatrix(REALSXP, maxHistLen, freeVars));
1603                         memcpy(REAL(Rdiff), diffWork.data(), sizeof(double) * maxHistLen * freeVars);
1604                         dbg.add("semDiff", Rdiff);
1605                 }
1606
1607                 if (paramHistLen.size()) {
1608                         SEXP Rphl;
1609                         Rf_protect(Rphl = Rf_allocVector(INTSXP, freeVars));
1610                         memcpy(INTEGER(Rphl), paramHistLen.data(), sizeof(int) * freeVars);
1611                         dbg.add("paramHistLen", Rphl);
1612                 }
1613
1614                 if (inputInfoMatrix) dbg.add("inputInfo", inputInfoMatrix);
1615                 if (outputInfoMatrix) dbg.add("outputInfo", outputInfoMatrix);
1616                 if (rateMatrix) dbg.add("rateMatrix", rateMatrix);
1617                 if (origEigenvalues) dbg.add("origEigenvalues", origEigenvalues);
1618
1619                 slots->add("debug", dbg.asR());
1620         }
1621 }
1622
1623 ComputeEM::~ComputeEM()
1624 {
1625         for (size_t rx=0; rx < ramsay.size(); ++rx) {
1626                 delete ramsay[rx];
1627         }
1628         ramsay.clear();
1629
1630         delete fit1;
1631
1632         for (size_t hx=0; hx < estHistory.size(); ++hx) {
1633                 delete [] estHistory[hx];
1634         }
1635         estHistory.clear();
1636 }
1637
1638 enum ComputeInfoMethod omxCompute::stringToInfoMethod(const char *iMethod)
1639 {
1640         enum ComputeInfoMethod infoMethod = INFO_METHOD_DEFAULT; // to avoid gcc warning
1641         if (strcmp(iMethod, "sandwich")==0) {
1642                 infoMethod = INFO_METHOD_SANDWICH;
1643         } else if (strcmp(iMethod, "meat")==0) {
1644                 infoMethod = INFO_METHOD_MEAT;
1645         } else if (strcmp(iMethod, "bread")==0) {
1646                 infoMethod = INFO_METHOD_BREAD;
1647         } else if (strcmp(iMethod, "hessian")==0) {
1648                 infoMethod = INFO_METHOD_HESSIAN;
1649         } else {
1650                 Rf_error("Unknown information matrix estimation method '%s'", iMethod);
1651         }
1652         return infoMethod;
1653 }
1654
1655 void omxComputeOnce::initFromFrontend(SEXP rObj)
1656 {
1657         super::initFromFrontend(rObj);
1658
1659         SEXP slotValue;
1660         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("from")));
1661         for (int wx=0; wx < Rf_length(slotValue); ++wx) {
1662                 int objNum = INTEGER(slotValue)[wx];
1663                 if (objNum >= 0) {
1664                         omxMatrix *algebra = globalState->algebraList[objNum];
1665                         if (algebra->fitFunction) {
1666                                 setFreeVarGroup(algebra->fitFunction, varGroup);
1667                                 omxCompleteFitFunction(algebra);
1668                         }
1669                         algebras.push_back(algebra);
1670                 } else {
1671                         omxExpectation *expectation = globalState->expectationList[~objNum];
1672                         setFreeVarGroup(expectation, varGroup);
1673                         omxCompleteExpectation(expectation);
1674                         expectations.push_back(expectation);
1675                 }
1676         }
1677
1678         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("verbose")));
1679         verbose = Rf_asInteger(slotValue);
1680
1681         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("what")));
1682         int whatLen = Rf_length(slotValue);
1683         if (algebras.size()) {
1684                 for (int wx=0; wx < whatLen; ++wx) {
1685                         SEXP elem;
1686                         Rf_protect(elem = STRING_ELT(slotValue, wx));
1687                         const char *what = CHAR(elem);
1688                         if      (strcmp(what, "maxAbsChange")==0) mac = true;
1689                         else if (strcmp(what, "set-starting-values")==0) starting = true;
1690                         else if (strcmp(what, "fit")         ==0) fit = true;
1691                         else if (strcmp(what, "gradient")    ==0) gradient = true;
1692                         else if (strcmp(what, "hessian")     ==0) hessian = true;
1693                         else if (strcmp(what, "information") ==0) infoMat = true;
1694                         else if (strcmp(what, "ihessian")    ==0) ihessian = true;
1695                         else omxRaiseErrorf("mxComputeOnce: don't know how to compute %s", what);
1696                 }
1697
1698                 if (hessian && infoMat) Rf_error("Cannot compute the Hessian and Fisher Information matrix simultaneously");
1699         } else {
1700                 for (int wx=0; wx < whatLen; ++wx) {
1701                         SEXP elem;
1702                         Rf_protect(elem = STRING_ELT(slotValue, wx));
1703                         predict.push_back(CHAR(elem));
1704                 }
1705         }
1706
1707         Rf_protect(slotValue = R_do_slot(rObj, Rf_install(".is.bestfit")));
1708         isBestFit = Rf_asLogical(slotValue);
1709
1710         bool howConflict = false;
1711         Rf_protect(slotValue = R_do_slot(rObj, Rf_install("how")));
1712         if (Rf_length(slotValue) > 1) {
1713                 omxRaiseErrorf("mxComputeOnce: more than one method specified");
1714         } else if (Rf_length(slotValue) == 1) {
1715                 SEXP elem;
1716                 Rf_protect(elem = STRING_ELT(slotValue, 0));
1717                 if (algebras.size()) {
1718                         const char *iMethod = CHAR(elem);
1719                         if (infoMat) {
1720                                 infoMethod = stringToInfoMethod(iMethod);
1721                                 if (infoMethod == INFO_METHOD_MEAT && gradient && whatLen == 2) {
1722                                         //OK
1723                                 } else if (whatLen > 1) {
1724                                         howConflict = true;
1725                                 }
1726                         } else {
1727                                 omxRaiseErrorf("mxComputeOnce: unknown method %s requested", iMethod);
1728                         }
1729                 } else {
1730                         how = CHAR(elem);
1731                         if (whatLen > 1) howConflict = true;
1732                 }
1733         }
1734         if (howConflict) {
1735                 omxRaiseErrorf("mxComputeOnce: when how is specified, you can only compute one thing at a time");
1736         }
1737
1738         if (algebras.size() == 1 && algebras[0]->fitFunction) {
1739                 omxFitFunction *ff = algebras[0]->fitFunction;
1740                 if (gradient && !ff->gradientAvailable) {
1741                         Rf_error("Gradient requested but not available");
1742                 }
1743                 if ((hessian || ihessian || hgprod) && !ff->hessianAvailable) {
1744                         // add a separate flag for hgprod TODO
1745                         Rf_error("Hessian requested but not available");
1746                 }
1747                 // add check for information TODO
1748         }
1749 }
1750
1751 void omxComputeOnce::computeImpl(FitContext *fc)
1752 {
1753         if (algebras.size()) {
1754                 int want = 0;
1755                 if (starting) {
1756                         want |= FF_COMPUTE_STARTING;
1757                 }
1758                 if (mac) {
1759                         want |= FF_COMPUTE_MAXABSCHANGE;
1760                         fc->mac = 0;
1761                 }
1762                 if (fit) {
1763                         want |= FF_COMPUTE_FIT;
1764                         if (isBestFit) want |= FF_COMPUTE_BESTFIT;
1765                         fc->fit = 0;
1766                 }
1767                 if (gradient) {
1768                         want |= FF_COMPUTE_GRADIENT;
1769                         fc->grad = Eigen::VectorXd::Zero(fc->numParam);
1770                 }
1771                 if (hessian) {
1772                         want |= FF_COMPUTE_HESSIAN;
1773                         fc->clearHessian();
1774                 }
1775                 if (infoMat) {
1776                         want |= FF_COMPUTE_INFO;
1777                         fc->infoMethod = infoMethod;
1778                         fc->grad = Eigen::VectorXd::Zero(fc->numParam);
1779                         fc->clearHessian();
1780                         fc->preInfo();
1781                 }
1782                 if (ihessian) {
1783                         want |= FF_COMPUTE_IHESSIAN;
1784                         fc->clearHessian();
1785                 }
1786                 if (!want) return;
1787
1788                 for (size_t wx=0; wx < algebras.size(); ++wx) {
1789                         omxMatrix *algebra = algebras[wx];
1790                         if (algebra->fitFunction) {
1791                                 omxFitFunctionCompute(algebra->fitFunction, FF_COMPUTE_PREOPTIMIZE, fc);
1792                                 omxFitFunctionCompute(algebra->fitFunction, want, fc);
1793                                 fc->fit = algebra->data[0];
1794                                 if (infoMat) {
1795                                         fc->postInfo();
1796                                 }
1797                         } else {
1798                                 omxForceCompute(algebra);
1799                         }
1800                 }
1801         } else if (expectations.size()) {
1802                 if (predict.size() > 1) Rf_error("Not implemented");
1803                 for (size_t wx=0; wx < expectations.size(); ++wx) {
1804                         omxExpectation *expectation = expectations[wx];
1805                         omxExpectationCompute(expectation, predict[0], how);
1806                 }
1807         }
1808 }
1809
1810 void omxComputeOnce::reportResults(FitContext *fc, MxRList *slots, MxRList *out)
1811 {
1812         if (algebras.size()==0 || algebras[0]->fitFunction == NULL) return;
1813
1814         omxMatrix *algebra = algebras[0];
1815         omxPopulateFitFunction(algebra, out);
1816 }
1817
1818 void ComputeStandardError::reportResults(FitContext *fc, MxRList *slots, MxRList *)
1819 {
1820         fc->allocStderrs();  // at least report NAs
1821
1822         const size_t numParams = fc->numParam;
1823
1824         Eigen::VectorXd ihessDiag(fc->ihessDiag());
1825
1826         const double scale = fabs(Global->llScale);
1827
1828         // This function calculates the standard errors from the Hessian matrix
1829         // sqrt(scale * diag(solve(hessian)))
1830
1831         for(size_t i = 0; i < numParams; i++) {
1832                 double got = ihessDiag[i];
1833                 if (got <= 0) continue;
1834                 fc->stderrs[i] = sqrt(scale * got);
1835         }
1836 }
1837
1838 /*
1839 Date: Fri, 3 Jan 2014 14:02:34 -0600
1840 From: Michael Hunter <mhunter@ou.edu>
1841
1842 Determining positive definiteness of matrix is typically done by
1843 trying the Cholesky decomposition.  If it fails, the matrix is not
1844 positive definite; if it passes, the matrix is.  The benefit of the
1845 Cholesky is that it's much faster and easier to compute than a set of
1846 eigenvalues.
1847
1848 The BLAS/LAPACK routine DTRCO quickly computes a good approximation to the
1849 reciprocal condition number of a triangular matrix.  Hand it the Cholesky
1850 (a triangular matrix) the rest is history.  I don't think we need the
1851 exact condition number as long as it's just for finding very
1852 ill-conditioned problems.  For the solution to a linear system of
1853 equations, if you really care about the difference in precision between
1854 1e-14 and 1e-11, then the exact condition number is needed.  Otherwise, the
1855 approximation is faster and equally useful.
1856 */
1857 void ComputeHessianQuality::reportResults(FitContext *fc, MxRList *slots, MxRList *)
1858 {
1859         // See Luenberger & Ye (2008) Second Order Test (p. 190) and Condition Number (p. 239)
1860
1861         // memcmp is required here because NaN != NaN always
1862         if (fc->infoDefinite != NA_LOGICAL ||
1863             memcmp(&fc->infoCondNum, &NA_REAL, sizeof(double)) != 0) return; // already set elsewhere
1864
1865         int numParams = int(fc->varGroup->vars.size());
1866
1867         double *mat = fc->getDenseHessianish();
1868         omxBuffer<double> hessWork(numParams * numParams);
1869         memcpy(hessWork.data(), mat, sizeof(double) * numParams * numParams);
1870
1871         char jobz = 'N';
1872         char range = 'A';
1873         char uplo = 'U';
1874         double abstol = 0;
1875         int m;
1876         omxBuffer<double> w(numParams);
1877         double optWork;
1878         int lwork = -1;
1879         omxBuffer<int> iwork(5 * numParams);
1880         int info;
1881         double realIgn = 0;
1882         int intIgn = 0;
1883         F77_CALL(dsyevx)(&jobz, &range, &uplo, &numParams, hessWork.data(),
1884                          &numParams, &realIgn, &realIgn, &intIgn, &intIgn, &abstol, &m, w.data(),
1885                          NULL, &numParams, &optWork, &lwork, iwork.data(), NULL, &info);
1886
1887         lwork = optWork;
1888         omxBuffer<double> work(lwork);
1889         F77_CALL(dsyevx)(&jobz, &range, &uplo, &numParams, hessWork.data(),
1890                          &numParams, &realIgn, &realIgn, &intIgn, &intIgn, &abstol, &m, w.data(),
1891                          NULL, &numParams, work.data(), &lwork, iwork.data(), NULL, &info);
1892         if (info < 0) {
1893                 Rf_error("dsyevx %d", info);
1894         } else if (info) {
1895                 return;
1896         }
1897
1898         bool definite = true;
1899         bool neg = w[0] < 0;
1900         for (int px=1; px < numParams; ++px) {
1901                 if ((w[px] < 0) ^ neg) {
1902                         definite = false;
1903                         break;
1904                 }
1905         }
1906
1907         fc->infoDefinite = definite;
1908
1909         if (definite) {
1910                 double ev[2] = { fabs(w[0]), fabs(w[numParams-1]) };
1911                 if (ev[0] < ev[1]) std::swap(ev[0], ev[1]);
1912                 double got = ev[0] / ev[1];
1913                 if (std::isfinite(got)) fc->infoCondNum = got;
1914         }
1915 }
1916
1917 void ComputeReportDeriv::reportResults(FitContext *fc, MxRList *, MxRList *result)
1918 {
1919         size_t numFree = fc->numParam;
1920
1921         if (fc->wanted & FF_COMPUTE_GRADIENT) {
1922                 if (!fc->grad.data()) {
1923                         Rf_warning("ComputeReportDeriv: Gradient requested but not available");
1924                 } else {
1925                         SEXP Rgradient;
1926                         Rf_protect(Rgradient = Rf_allocVector(REALSXP, numFree));
1927                         memcpy(REAL(Rgradient), fc->grad.data(), sizeof(double) * numFree);
1928                         result->add("gradient", Rgradient);
1929                 }
1930         }
1931         if (fc->wanted & FF_COMPUTE_HESSIAN) {
1932                 SEXP Rhessian;
1933                 Rf_protect(Rhessian = Rf_allocMatrix(REALSXP, numFree, numFree));
1934                 fc->copyDenseHess(REAL(Rhessian));
1935                 result->add("hessian", Rhessian);
1936         }
1937         if (fc->wanted & FF_COMPUTE_IHESSIAN) {
1938                 SEXP Rihessian;
1939                 Rf_protect(Rihessian = Rf_allocMatrix(REALSXP, numFree, numFree));
1940                 fc->copyDenseIHess(REAL(Rihessian));
1941                 result->add("ihessian", Rihessian);
1942         }
1943 }