Reimagine CSOLNPFit as GradientOptimizerAPI
[openmx:openmx.git] / src / Compute.h
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 #ifndef _OMX_COMPUTE_H_
18 #define _OMX_COMPUTE_H_
19
20 #define R_NO_REMAP
21 #include <R.h>
22 #include <Rinternals.h>
23 #include "omxDefines.h"
24 #include "Eigen/SparseCore"
25 #include "glue.h"
26
27 // See R/MxRunHelperFunctions.R optimizerMessages
28 // Also see the NPSOL manual, section 6 "Error Indicators and Warnings"
29 // These are ordered from good to bad so we can use max() on a set
30 // of inform results to obtain a bound on convergence status.
31 typedef int ComputeInform;
32 #define INFORM_UNINITIALIZED -1
33 #define INFORM_CONVERGED_OPTIMUM 0
34 #define INFORM_UNCONVERGED_OPTIMUM 1
35         // The final iterate satisfies the optimality conditions to the accuracy requested,
36         // but the sequence of iterates has not yet converged.
37         // Optimizer was terminated because no further improvement
38         // could be made in the merit function (Mx status GREEN).
39 #define INFORM_LINEAR_CONSTRAINTS_INFEASIBLE 2
40         // The linear constraints and bounds could not be satisfied.
41         // The problem has no feasible solution.
42 #define INFORM_NONLINEAR_CONSTRAINTS_INFEASIBLE 3
43         // The nonlinear constraints and bounds could not be satisfied.
44         // The problem may have no feasible solution.
45 #define INFORM_ITERATION_LIMIT 4
46         // The major iteration limit was reached (Mx status BLUE).
47 #define INFORM_NOT_AT_OPTIMUM 6
48         // The model does not satisfy the first-order optimality conditions (i.e. gradient close to zero)
49         // to the required accuracy, and no improved point for the
50         // merit function could be found during the final linesearch (Mx status RED)
51 #define INFORM_BAD_DERIVATIVES 7
52         // The function derivates returned by funcon or funobj
53         // appear to be incorrect.
54 #define INFORM_INVALID_PARAM 9
55         // An input parameter was invalid'
56 #define INFORM_STARTING_VALUES_INFEASIBLE 10
57
58 enum ComputeInfoMethod {
59         INFO_METHOD_DEFAULT,
60         INFO_METHOD_HESSIAN,
61         INFO_METHOD_SANDWICH,
62         INFO_METHOD_BREAD,
63         INFO_METHOD_MEAT
64 };
65
66 struct HessianBlock {
67         std::vector<int> vars;  // global freeVar ID in order
68         Eigen::MatrixXd mat;    // vars * vars, only upper triangle referenced
69         Eigen::MatrixXd mmat;   // including subblocks
70         Eigen::MatrixXd imat;
71         std::vector< HessianBlock* > subBlocks;
72         bool merge;
73         int useId;
74
75         HessianBlock() : merge(false), useId(0) {}
76         HessianBlock *clone();
77         void addSubBlocks();
78         int estNonZero() const;
79 };
80
81 // The idea of FitContext is to eventually enable fitting from
82 // multiple starting values in parallel.
83
84 class FitContext {
85         static omxFitFunction *RFitFunction;
86
87         FitContext *parent;
88
89         std::vector<HessianBlock*> allBlocks;
90         std::vector<HessianBlock*> mergeBlocks;
91         std::vector<HessianBlock*> blockByVar;
92
93         bool haveSparseHess;
94         Eigen::SparseMatrix<double> sparseHess;
95         bool haveSparseIHess;
96         Eigen::SparseMatrix<double> sparseIHess;
97         int estNonZero;
98         int minBlockSize;
99         int maxBlockSize;
100
101         bool haveDenseHess;
102         Eigen::MatrixXd hess;
103         bool haveDenseIHess;
104         Eigen::MatrixXd ihess;
105
106         void init();
107         void analyzeHessian();
108         void analyzeHessianBlock(HessianBlock *hb);
109         void testMerge();
110
111         std::string IterationError;
112
113  public:
114         FreeVarGroup *varGroup;
115         omxState *state;
116         size_t numParam;               // change to int type TODO
117         std::vector<int> mapToParent;
118         double mac;
119         double fit;
120         double *est;
121         std::vector<const char*> flavor;
122         Eigen::VectorXd grad;
123         int infoDefinite;
124         double infoCondNum;
125         double *stderrs;   // plural to distinguish from stdio's stderr
126         enum ComputeInfoMethod infoMethod;
127         double *infoA; // sandwich, the bread
128         double *infoB; // sandwich, the meat
129         int iterations;
130         ComputeInform inform;
131         int wanted;
132         std::vector< class FitContext* > childList;
133
134         FitContext(omxState *_state, std::vector<double> &startingValues);
135         FitContext(FitContext *parent, FreeVarGroup *group);
136         void createChildren();
137         void allocStderrs();
138         void copyParamToModel();
139         void copyParamToModelClean();
140         double *take(int want);
141         omxMatrix *lookupDuplicate(omxMatrix* element);
142         void maybeCopyParamToModel(omxState* os);
143         void updateParent();
144         void updateParentAndFree();
145         template <typename T> void moveInsideBounds(std::vector<T> &prevEst);
146         void log(int what);
147         ~FitContext();
148         
149         // deriv related
150         void clearHessian();
151         void negateHessian();
152         void queue(HessianBlock *hb);
153         void refreshDenseHess();
154         void refreshDenseIHess();
155         void refreshSparseHess();
156         bool refreshSparseIHess(); // NOTE: produces an ihess with filtered eigenvalues
157         Eigen::VectorXd ihessGradProd();
158         double *getDenseHessUninitialized();
159         double *getDenseIHessUninitialized();
160         double *getDenseHessianish();  // either a Hessian or inverse Hessian, remove TODO
161         void copyDenseHess(double *dest);
162         void copyDenseIHess(double *dest);
163         Eigen::VectorXd ihessDiag();
164         void preInfo();
165         void postInfo();
166         void resetIterationError();
167         void recordIterationError(const char* msg, ...) __attribute__((format (printf, 2, 3)));
168
169         std::string getIterationError();
170
171         static void cacheFreeVarDependencies();
172         static void setRFitFunction(omxFitFunction *rff);
173 };
174
175 typedef std::vector< std::pair<int, MxRList*> > LocalComputeResult;
176
177 class omxCompute {
178         int computeId;
179  protected:
180         virtual void reportResults(FitContext *fc, MxRList *slots, MxRList *glob) {};
181         void collectResultsHelper(FitContext *fc, std::vector< omxCompute* > &clist,
182                                   LocalComputeResult *lcr, MxRList *out);
183         static enum ComputeInfoMethod stringToInfoMethod(const char *iMethod);
184  public:
185         const char *name;
186         FreeVarGroup *varGroup;
187         omxCompute();
188         virtual void initFromFrontend(omxState *, SEXP rObj);
189         void compute(FitContext *fc);
190         virtual void computeImpl(FitContext *fc) {}
191         virtual void collectResults(FitContext *fc, LocalComputeResult *lcr, MxRList *out);
192         virtual ~omxCompute();
193 };
194
195 omxCompute *omxNewCompute(omxState* os, const char *type);
196
197 omxCompute *newComputeGradientDescent();
198 omxCompute *newComputeNumericDeriv();
199 omxCompute *newComputeNewtonRaphson();
200 omxCompute *newComputeConfidenceInterval();
201
202 void omxApproxInvertPosDefTriangular(int dim, double *hess, double *ihess, double *stress);
203 void omxApproxInvertPackedPosDefTriangular(int dim, int *mask, double *packedHess, double *stress);
204 SEXP sparseInvert_wrapper(SEXP mat);
205
206 struct CSOLNPFit {
207         const char *optName;
208         FitContext *fc;
209
210         Eigen::VectorXd solLB;
211         Eigen::VectorXd solUB;
212
213         Eigen::VectorXd solIneqLB;
214         Eigen::VectorXd solIneqUB;
215
216         Eigen::VectorXd equality;
217         Eigen::VectorXd inequality;
218
219         // output
220         int informOut;
221         double fitOut;
222         Eigen::VectorXd gradOut;
223         Eigen::MatrixXd hessOut;
224
225         CSOLNPFit(const char *optName, FitContext *fc) : optName(optName), fc(fc) {}
226         virtual double solFun(double *myPars, int* mode, int verbose) = 0;
227         virtual void solEqBFun(int verbose) = 0;
228         virtual void myineqFun(int verbose) = 0;
229 };
230
231 struct RegularFit : CSOLNPFit {
232         typedef CSOLNPFit super;
233         omxMatrix *fitMatrix;
234         void setupIneqConstraintBounds();
235
236         RegularFit(const char *optName, FitContext *fc, omxMatrix *fmat);
237         virtual double solFun(double *myPars, int* mode, int verbose);
238         virtual void solEqBFun(int verbose);
239         virtual void myineqFun(int verbose);
240 };
241
242 struct ConfidenceIntervalFit : RegularFit {
243         typedef RegularFit super;
244         int currentInterval;
245         bool calcLower;
246
247         ConfidenceIntervalFit(const char *optName, FitContext *fc, omxMatrix *fmat, int curInt, bool lower);
248         virtual double solFun(double *myPars, int* mode, int verbose);
249 };
250
251 #endif