NLOPT and Simulated annealing added
[openmx:openmx.git] / src / ComputeGD.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 // Maybe integrate http://ab-initio.mit.edu/wiki/index.php/NLopt_Introduction ?
18
19 #include "omxDefines.h"
20 #include "omxState.h"
21 #include "omxFitFunction.h"
22 #include "omxNPSOLSpecific.h"
23 #include "omxExportBackendState.h"
24 #include "omxCsolnp.h"
25 #include "nloptcpp.h"
26 #include "Compute.h"
27 #include "npsolswitch.h"
28 #include "glue.h"
29
30 enum OptEngine {
31         OptEngine_NPSOL,
32         OptEngine_CSOLNP,
33     OptEngine_NLOPT
34 };
35
36 class ComputeGDBase : public omxCompute {
37 protected:
38         typedef omxCompute super;
39         enum OptEngine engine;
40         omxMatrix *fitMatrix;
41         int verbose;
42         double optimalityTolerance;
43
44         virtual void initFromFrontend(SEXP rObj);
45 };
46
47 class omxComputeGD : public ComputeGDBase {
48         typedef ComputeGDBase super;
49         bool useGradient;
50         SEXP hessChol;
51
52         int warmStartSize;
53         double *warmStart;
54     
55 public:
56         omxComputeGD();
57         virtual void initFromFrontend(SEXP rObj);
58         virtual void computeImpl(FitContext *fc);
59         virtual void reportResults(FitContext *fc, MxRList *slots, MxRList *out);
60 };
61
62 class omxCompute *newComputeGradientDescent()
63 {
64         return new omxComputeGD();
65 }
66
67 class ComputeCI : public ComputeGDBase {
68         typedef ComputeGDBase super;
69         SEXP intervals, intervalCodes;
70
71 public:
72         ComputeCI();
73         virtual void initFromFrontend(SEXP rObj);
74         virtual void computeImpl(FitContext *fc);
75         virtual void reportResults(FitContext *fc, MxRList *slots, MxRList *out);
76 };
77
78 omxCompute *newComputeConfidenceInterval()
79 {
80         return new ComputeCI();
81 }
82
83 omxComputeGD::omxComputeGD()
84 {
85         hessChol = NULL;
86         warmStart = NULL;
87 }
88
89 void ComputeGDBase::initFromFrontend(SEXP rObj)
90 {
91         super::initFromFrontend(rObj);
92
93         SEXP slotValue;
94         fitMatrix = omxNewMatrixFromSlot(rObj, globalState, "fitfunction");
95         setFreeVarGroup(fitMatrix->fitFunction, varGroup);
96         omxCompleteFitFunction(fitMatrix);
97
98         ScopedProtect p1(slotValue, R_do_slot(rObj, Rf_install("verbose")));
99         verbose = Rf_asInteger(slotValue);
100     
101         ScopedProtect p2(slotValue, R_do_slot(rObj, Rf_install("tolerance")));
102         optimalityTolerance = Rf_asReal(slotValue);
103     
104         ScopedProtect p3(slotValue, R_do_slot(rObj, Rf_install("engine")));
105         const char *engine_name = CHAR(Rf_asChar(slotValue));
106         if (strcmp(engine_name, "CSOLNP")==0) {
107                 engine = OptEngine_CSOLNP;
108         } else if (strcmp(engine_name, "NLOPT")==0) {
109         engine = OptEngine_NLOPT;
110     } else if (strcmp(engine_name, "NPSOL")==0) {
111 #if HAS_NPSOL
112                 engine = OptEngine_NPSOL;
113 #else
114                 Rf_error("NPSOL is not available in this build");
115 #endif
116         } else {
117                 Rf_error("%s: engine %s unknown", name, engine_name);
118         }
119 }
120
121 void omxComputeGD::initFromFrontend(SEXP rObj)
122 {
123         super::initFromFrontend(rObj);
124     
125         SEXP slotValue;
126         ScopedProtect p1(slotValue, R_do_slot(rObj, Rf_install("useGradient")));
127         if (Rf_length(slotValue)) {
128                 useGradient = Rf_asLogical(slotValue);
129         } else {
130                 useGradient = Global->analyticGradients;
131         }
132
133         ScopedProtect p2(slotValue, R_do_slot(rObj, Rf_install("warmStart")));
134         if (!Rf_isNull(slotValue)) {
135                 SEXP matrixDims;
136                 Rf_protect(matrixDims = Rf_getAttrib(slotValue, R_DimSymbol));
137                 int *dimList = INTEGER(matrixDims);
138                 int rows = dimList[0];
139                 int cols = dimList[1];
140                 if (rows != cols) Rf_error("%s: warmStart matrix must be square", name);
141
142                 warmStartSize = rows;
143                 warmStart = REAL(slotValue);
144         }
145 }
146
147 void omxComputeGD::computeImpl(FitContext *fc)
148 {
149     size_t numParam = varGroup->vars.size();
150         if (numParam <= 0) {
151                 omxRaiseErrorf("%s: model has no free parameters", name);
152                 return;
153         }
154     
155         omxFitFunctionCompute(fitMatrix->fitFunction, FF_COMPUTE_PREOPTIMIZE, fc);
156
157         fc->createChildren();
158     
159         switch (engine) {
160         case OptEngine_NPSOL:{
161 #if HAS_NPSOL
162                 if (!hessChol) {
163                         Rf_protect(hessChol = Rf_allocMatrix(REALSXP, numParam, numParam));
164                 }
165                 bool doWarm = false;
166                 if (warmStart) {
167                         if (warmStartSize != int(numParam)) {
168                                 Rf_warning("%s: warmStart size %d does not match number of free parameters %d (ignored)",
169                                            warmStartSize, numParam);
170                         } else {
171                                 memcpy(REAL(hessChol), warmStart, sizeof(double) * numParam * numParam);
172                                 doWarm = true;
173                         }
174                 }
175                 omxInvokeNPSOL(fitMatrix, fc, &fc->inform, useGradient, varGroup, verbose,
176                                REAL(hessChol), optimalityTolerance, doWarm);
177                 Eigen::Map<Eigen::MatrixXd> hc(REAL(hessChol), numParam, numParam);
178                 Eigen::MatrixXd hcT = hc.transpose();
179                 Eigen::Map<Eigen::MatrixXd> dest(fc->getDenseHessUninitialized(), numParam, numParam);
180                 dest.noalias() = hcT * hc;
181 #endif
182                 break;}
183         case OptEngine_CSOLNP:
184             omxInvokeCSOLNP(fitMatrix, fc, &fc->inform, varGroup, verbose,
185                             fc->getDenseHessUninitialized(), optimalityTolerance);
186         case OptEngine_NLOPT:
187             omxInvokeNLOPTorSANN(fitMatrix, fc, &fc->inform, varGroup, verbose,
188                 fc->getDenseHessUninitialized(), optimalityTolerance);
189             break;
190         default: Rf_error("huh?");
191         }
192         fc->wanted |= FF_COMPUTE_GRADIENT | FF_COMPUTE_HESSIAN;
193     
194         // It seems we cannot avoid this. Both optimizers can terminate
195         // with fit not at the optimum.
196         ComputeFit(fitMatrix, FF_COMPUTE_FIT, fc);
197
198         if (fitMatrix->rows == 1) {
199                 if (!std::isfinite(fc->fit) || fc->fit == 1e24) {  // remove magic number 1e24 TODO
200                         std::string diag = fc->getIterationError();
201                         omxRaiseErrorf("MxComputeGradientDescent: fitfunction %s evaluated to %f (%s)",
202                                        fitMatrix->name, fc->fit, diag.c_str());
203                         return;
204                 }
205         }
206
207         fc->wanted |= FF_COMPUTE_BESTFIT;
208     /*printf("fc->hess in computeGD\n");
209     printf("%2f", fc->hess[0]); putchar('\n');
210     printf("%2f", fc->hess[1]); putchar('\n');
211     printf("%2f", fc->hess[2]); putchar('\n');
212     */
213 }
214
215 void omxComputeGD::reportResults(FitContext *fc, MxRList *slots, MxRList *out)
216 {
217         omxPopulateFitFunction(fitMatrix, out);
218     
219     /*printf("fc->hess in computeGD:report results\n");
220     printf("%2f", fc->hess[0]); putchar('\n');
221     printf("%2f", fc->hess[1]); putchar('\n');
222     printf("%2f", fc->hess[2]); putchar('\n');
223 */
224         if (engine == OptEngine_NPSOL) {
225                 out->add("hessianCholesky", hessChol);
226         }
227 }
228
229 ComputeCI::ComputeCI()
230 {
231         intervals = 0;
232         intervalCodes = 0;
233 }
234
235 void ComputeCI::initFromFrontend(SEXP rObj)
236 {
237         super::initFromFrontend(rObj);
238 }
239
240 void ComputeCI::computeImpl(FitContext *fc)
241 {
242         Global->unpackConfidenceIntervals();
243
244         int numInts = (int) Global->intervalList.size();
245         if (verbose >= 1) mxLog("%s: starting work on %d intervals", name, numInts);
246         if (!numInts) return;
247
248         // I'm not sure why INFORM_NOT_AT_OPTIMUM is okay, but that's how it was.
249         if (fc->inform >= INFORM_LINEAR_CONSTRAINTS_INFEASIBLE && fc->inform != INFORM_NOT_AT_OPTIMUM) {
250                 // TODO: allow forcing
251                 Rf_warning("Not calculating confidence intervals because of optimizer status %d", fc->inform);
252                 return;
253         }
254
255         Eigen::ArrayXd mle(fc->numParam);
256         memcpy(mle.data(), fc->est, sizeof(double) * fc->numParam);
257
258         Rf_protect(intervals = Rf_allocMatrix(REALSXP, numInts, 3));
259         Rf_protect(intervalCodes = Rf_allocMatrix(INTSXP, numInts, 2));
260
261         switch (engine) {
262         case OptEngine_NPSOL:
263 #if HAS_NPSOL
264                 omxNPSOLConfidenceIntervals(fitMatrix, fc, optimalityTolerance);
265 #endif
266                 break;
267         case OptEngine_CSOLNP:
268                 omxCSOLNPConfidenceIntervals(fitMatrix, fc, verbose, optimalityTolerance);
269                 break;
270         default:
271                 Rf_error("huh?");
272         }
273
274         if(OMX_DEBUG) { mxLog("Populating CIs for %d fit functions.", numInts); }
275
276         memcpy(fc->est, mle.data(), sizeof(double) * fc->numParam);
277         fc->copyParamToModel();
278
279         Eigen::Map< Eigen::ArrayXXd > interval(REAL(intervals), numInts, 3);
280         int* intervalCode = INTEGER(intervalCodes);
281         for(int j = 0; j < numInts; j++) {
282                 omxConfidenceInterval *oCI = Global->intervalList[j];
283                 interval(j, 0) = oCI->min;
284                 omxRecompute(oCI->matrix, FF_COMPUTE_FIT, fc);
285                 interval(j, 1) = omxMatrixElement(oCI->matrix, oCI->row, oCI->col);
286                 interval(j, 2) = oCI->max;
287                 intervalCode[j] = oCI->lCode;
288                 intervalCode[j + numInts] = oCI->uCode;
289         }
290 }
291
292 void ComputeCI::reportResults(FitContext *fc, MxRList *slots, MxRList *out)
293 {
294         if (!intervals) return;
295
296         int numInt = (int) Global->intervalList.size();
297
298         SEXP dimnames;
299         SEXP names;
300         Rf_protect(dimnames = Rf_allocVector(VECSXP, 2));
301         Rf_protect(names = Rf_allocVector(STRSXP, 3));
302         SET_STRING_ELT(names, 0, Rf_mkChar("lbound"));
303         SET_STRING_ELT(names, 1, Rf_mkChar("estimate"));
304         SET_STRING_ELT(names, 2, Rf_mkChar("ubound"));
305         SET_VECTOR_ELT(dimnames, 1, names);
306
307         Rf_protect(names = Rf_allocVector(STRSXP, numInt)); //shared between the two matrices
308         for (int nx=0; nx < numInt; ++nx) {
309                 omxConfidenceInterval *ci = Global->intervalList[nx];
310                 SET_STRING_ELT(names, nx, Rf_mkChar(ci->name));
311         }
312         SET_VECTOR_ELT(dimnames, 0, names);
313
314         Rf_setAttrib(intervals, R_DimNamesSymbol, dimnames);
315
316         out->add("confidenceIntervals", intervals);
317
318         Rf_protect(dimnames = Rf_allocVector(VECSXP, 2));
319         SET_VECTOR_ELT(dimnames, 0, names);
320
321         Rf_protect(names = Rf_allocVector(STRSXP, 2));
322         SET_STRING_ELT(names, 0, Rf_mkChar("lbound"));
323         SET_STRING_ELT(names, 1, Rf_mkChar("ubound"));
324         SET_VECTOR_ELT(dimnames, 1, names);
325
326         Rf_setAttrib(intervalCodes, R_DimNamesSymbol, dimnames);
327
328         out->add("confidenceIntervalCodes", intervalCodes);
329 }