Reimagine CSOLNPFit as GradientOptimizerAPI
[openmx:openmx.git] / src / omxCsolnp.cpp
1 /*
2  *  Copyright 2007-2012 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 <ctype.h>
18 #include <limits>
19 #define R_NO_REMAP
20 #include <R.h>
21 #include <Rinternals.h>
22 #include "omxState.h"
23 #include "omxNPSOLSpecific.h"
24 #include "omxMatrix.h"
25 #include "glue.h"
26 #include "omxImportFrontendState.h"
27 #include "omxCsolnp.h"
28 #include "omxBuffer.h"
29
30 static int majIter = 400;
31 static int minIter = 800;
32 static double funcPrecision = 1.0e-7;
33
34 void omxInvokeCSOLNP(omxMatrix *fitMatrix, FitContext *fc,
35                      int *inform_out, FreeVarGroup *freeVarGroup,
36                      int verbose, double *hessOut, double tolerance)
37 {
38     fc->grad.resize(fc->numParam);
39     
40     int n = int(freeVarGroup->vars.size());
41     
42     Eigen::Map< Eigen::VectorXd > myPars(fc->est, n);
43     
44     Eigen::Array<double, 5, 1> myControl;
45     myControl[0] = 1.0;
46     myControl[1] = majIter;
47     myControl[2] = minIter;
48     myControl[3] = funcPrecision;
49     myControl[4] = std::isfinite(tolerance)? tolerance : 1.0e-9;
50     
51     RegularFit rf("CSOLNP", fc, fitMatrix);
52     solnp(myPars.data(), rf, myControl, verbose);
53     
54     fc->fit = rf.fitOut;
55     if (verbose >= 1) {
56         mxLog("final objective value is: \n");
57         mxLog("%2f", fc->fit);
58     }
59     
60     if (verbose>= 1){
61         mxLog("final myPars value is: \n");
62         for (int i = 0; i < myPars.size(); i++) mxLog("%f", myPars[i]);
63     }
64     
65     *inform_out = rf.informOut;
66     
67     if (rf.gradOut.size()) {
68             fc->grad = rf.gradOut.tail(n);
69             Eigen::Map< Eigen::MatrixXd > hess(hessOut, n, n);
70             hess = rf.hessOut.bottomRightCorner(n, n);
71     }
72     
73     fc->copyParamToModel();
74 }
75
76 // Mostly duplicated code in omxNPSOLConfidenceIntervals
77 // needs to be refactored so there is only 1 copy of CI
78 // code that can use whatever optimizer is provided.
79 void omxCSOLNPConfidenceIntervals(omxMatrix *fitMatrix, FitContext *opt, int verbose, double tolerance)
80 {
81     const int ciMaxIterations = Global->ciMaxIterations;
82     FitContext fc(opt, opt->varGroup);
83     fc.createChildren();
84     FreeVarGroup *freeVarGroup = fc.varGroup;
85     
86     const int n = int(freeVarGroup->vars.size());
87     
88     Eigen::Array<double, 5, 1> myControl;
89     myControl[0] = 1.0;
90     myControl[1] = majIter;
91     myControl[2] = minIter;
92     myControl[3] = funcPrecision;
93     myControl[4] = std::isfinite(tolerance)? tolerance : 1.0e-16;
94
95     if(OMX_DEBUG) { mxLog("Calculating likelihood-based confidence intervals."); }
96     
97     const double objDiff = 1.e-4;     // TODO : Use function precision to determine CI jitter?
98     
99     for(int i = 0; i < (int) Global->intervalList.size(); i++) {
100         omxConfidenceInterval *currentCI = Global->intervalList[i];
101         
102         const char *matName = "anonymous matrix";
103         if (currentCI->matrix->name) {
104             matName = currentCI->matrix->name;
105         }
106         
107         currentCI->lbound += opt->fit;          // Convert from offsets to targets
108         currentCI->ubound += opt->fit;          // Convert from offsets to targets
109         
110         for (int lower=0; lower <= 1; ++lower) {
111                 if (lower  && !std::isfinite(currentCI->lbound)) continue;
112                 if (!lower && !std::isfinite(currentCI->ubound)) continue;
113
114                 memcpy(fc.est, opt->est, n * sizeof(double)); // Reset to previous optimum
115         
116                 int tries = 0;
117                 int inform = -1;
118                 double bestFit = std::numeric_limits<double>::max();
119             
120                 while (inform!= 0 && ++tries <= ciMaxIterations) {
121                         Global->checkpointMessage(opt, opt->est, "%s[%d, %d] %s CI (try %d)",
122                                                   matName, currentCI->row + 1, currentCI->col + 1,
123                                                   lower? "lower" : "upper", tries);
124
125                         ConfidenceIntervalFit cif("CSOLNP", &fc, fitMatrix, i, lower);
126                         solnp(fc.est, cif, myControl, verbose);
127                 
128                         if(cif.fitOut < bestFit) {
129                                 double val = omxMatrixElement(currentCI->matrix, currentCI->row, currentCI->col);
130                                 if (lower) currentCI->min = val;
131                                 else       currentCI->max = val;
132                                 bestFit = cif.fitOut;
133                         }
134
135                         inform = cif.informOut;
136                         if (lower) currentCI->lCode = inform;
137                         else       currentCI->uCode = inform;
138                         if(verbose>=1) { mxLog("inform(%d,%d) is: %d", i, lower, inform);}
139                         if(inform == 0) break;
140
141                         bool jitter = TRUE;
142                         for(int j = 0; j < n; j++) {
143                                 if(fabs(fc.est[j] - opt->est[j]) > objDiff) {
144                                         jitter = FALSE;
145                                         break;
146                                 }
147                         }
148                         if(jitter) {
149                                 for(int j = 0; j < n; j++) {
150                                         double sign = 2 * (tries % 2) - 1;
151                                         fc.est[j] = opt->est[j] + sign * objDiff * tries;
152                                 }
153                         }
154                 }
155         }
156     }
157 }
158
159 void CSOLNPOpt_majIter(const char *optionValue)
160 {
161     majIter = atoi(optionValue);
162 }
163
164 void CSOLNPOpt_minIter(const char *optionValue)
165 {
166     minIter = atoi(optionValue);
167 }
168
169 void CSOLNPOpt_FuncPrecision(const char *optionValue)
170 {
171     funcPrecision = atof(optionValue);
172 }