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