Remove most instances of setFinalReturns
[openmx:openmx.git] / src / ComputeNR.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 "omxState.h"
18 #include "omxFitFunction.h"
19 #include "omxExportBackendState.h"
20 #include "Compute.h"
21
22 class ComputeNR : public omxComputeOperation {
23         typedef omxComputeOperation super;
24         omxMatrix *fitMatrix;
25
26         int maxIter;
27         double tolerance;
28         int inform, iter;
29
30 public:
31         ComputeNR();
32         virtual void initFromFrontend(SEXP rObj);
33         virtual void compute(FitContext *fc);
34         virtual void reportResults(FitContext *fc, MxRList *out);
35         virtual double getOptimizerStatus() { return inform; }  // backward compatibility
36 };
37
38 class omxCompute *newComputeNewtonRaphson()
39 {
40         return new ComputeNR();
41 }
42
43 ComputeNR::ComputeNR()
44 {
45         inform = 0;
46         iter = 0;
47 }
48
49 void ComputeNR::initFromFrontend(SEXP rObj)
50 {
51         super::initFromFrontend(rObj);
52
53         fitMatrix = omxNewMatrixFromSlot(rObj, globalState, "fitfunction");
54         setFreeVarGroup(fitMatrix->fitFunction, varGroup);
55         omxCompleteFitFunction(fitMatrix);
56
57         if (!fitMatrix->fitFunction->hessianAvailable ||
58             !fitMatrix->fitFunction->gradientAvailable) {
59                 error("Newton-Raphson requires derivatives");
60         }
61
62         SEXP slotValue;
63         PROTECT(slotValue = GET_SLOT(rObj, install("maxIter")));
64         maxIter = INTEGER(slotValue)[0];
65
66         PROTECT(slotValue = GET_SLOT(rObj, install("tolerance")));
67         tolerance = REAL(slotValue)[0];
68         if (tolerance <= 0) error("tolerance must be positive");
69 }
70
71 void ComputeNR::compute(FitContext *fc)
72 {
73         // complain if there are non-linear constraints TODO
74
75         size_t numParam = varGroup->vars.size();
76         if (numParam <= 0) {
77                 error("Model has no free parameters");
78                 return;
79         }
80
81         omxFitFunctionCompute(fitMatrix->fitFunction, FF_COMPUTE_PREOPTIMIZE, fc);
82
83         iter = 0;
84         double prevLL = nan("unset");
85         bool decreasing = TRUE;
86
87         while (1) {
88                 const int want = FF_COMPUTE_FIT|FF_COMPUTE_GRADIENT|FF_COMPUTE_HESSIAN;
89
90                 OMXZERO(fc->grad, numParam);
91                 OMXZERO(fc->hess, numParam * numParam);
92
93                 omxFitFunctionCompute(fitMatrix->fitFunction, want, fc);
94
95                 // Only need LL for diagnostics; Can avoid computing it? TODO
96                 double LL = fitMatrix->data[0];
97                 if (isfinite(prevLL) && prevLL < LL - tolerance) decreasing = FALSE;
98                 prevLL = LL;
99
100                 fc->fixHessianSymmetry();
101                 //              fc->log(FF_COMPUTE_ESTIMATE|FF_COMPUTE_GRADIENT|FF_COMPUTE_HESSIAN);
102
103                 std::vector<double> ihess(numParam * numParam);
104                 for (size_t rx=0; rx < numParam; rx++) {
105                         for (size_t cx=0; cx < numParam; cx++) {
106                                 ihess[rx * numParam + cx] = rx==cx? 1 : 0;
107                         }
108                 }
109                 std::vector<int> ipiv(numParam);
110                 int info;
111                 int dim = int(numParam);
112                 F77_CALL(dgesv)(&dim, &dim, fc->hess, &dim, ipiv.data(),
113                                 ihess.data(), &dim, &info);
114                 if (info < 0) error("Arg %d is invalid", -info);
115                 if (info > 0) {
116                         omxRaiseErrorf(globalState, "Hessian is singular and cannot be inverted");
117                         break;
118                 }
119
120                 std::vector<double> adj(numParam);
121                 char trans = 'N';
122                 double alpha = -1;
123                 int incx = 1;
124                 double beta = 0;
125                 F77_CALL(dgemv)(&trans, &dim, &dim, &alpha, ihess.data(), &dim,
126                                 fc->grad, &incx, &beta, adj.data(), &incx);
127
128                 double maxAdj = 0;
129                 for (size_t px=0; px < numParam; ++px) {
130                         double param = fc->est[px];
131                         param += adj[px];
132                         omxFreeVar *fv = fc->varGroup->vars[px];
133                         if (param < fv->lbound) param = fv->lbound;
134                         if (param > fv->ubound) param = fv->ubound;
135                         double adj = fabs(param - fc->est[px]);
136                         if (maxAdj < adj)
137                                 maxAdj = adj;
138                         fc->est[px] = param;
139                 }
140                 fc->copyParamToModel(globalState);
141                 R_CheckUserInterrupt();
142                 if (maxAdj < tolerance || ++iter > maxIter) break;
143         }
144
145         // The check is too dependent on numerical precision to enable by default.
146         // Anyway, it's just a tool for developers.
147         if (0 && !decreasing) warning("Newton-Raphson iterations did not converge");
148 }
149
150 void ComputeNR::reportResults(FitContext *fc, MxRList *out)
151 {
152         if (Global->numIntervals) {
153                 warning("Confidence intervals are not implemented for Newton-Raphson");
154         }  
155
156         omxPopulateFitFunction(fitMatrix, out);
157
158         size_t numFree = varGroup->vars.size();
159
160         SEXP estimate;
161         PROTECT(estimate = allocVector(REALSXP, numFree));
162         memcpy(REAL(estimate), fc->est, sizeof(double) * numFree);
163
164         out->push_back(std::make_pair(mkChar("minimum"), ScalarReal(fc->fit)));
165         out->push_back(std::make_pair(mkChar("Minus2LogLikelihood"), ScalarReal(fc->fit)));
166         out->push_back(std::make_pair(mkChar("estimate"), estimate));
167
168         SEXP code, iterations;
169
170         PROTECT(code = NEW_NUMERIC(1));
171         REAL(code)[0] = inform;
172         out->push_back(std::make_pair(mkChar("nr.code"), code));
173
174         PROTECT(iterations = NEW_NUMERIC(1));
175         REAL(iterations)[0] = iter;
176         out->push_back(std::make_pair(mkChar("nr.iterations"), iterations));
177 }