Generic ComputeEM with Ramsay (1975) acceleration
[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 #include <R.h>
21 #include <Rinternals.h>
22
23 #include "types.h"
24
25 // See R/MxRunHelperFunctions.R npsolMessages
26 enum ComputeInform {
27         INFORM_UNINITIALIZED = -1,
28         INFORM_CONVERGED_OPTIMUM = 0,
29         INFORM_UNCONVERGED_OPTIMUM = 1,
30         // The final iterate satisfies the optimality conditions to the accuracy requested,
31         // but the sequence of iterates has not yet converged.
32         // Optimizer was terminated because no further improvement
33         // could be made in the merit function (Mx status GREEN).
34         INFORM_LINEAR_CONSTRAINTS_INFEASIBLE = 2,
35         // The linear constraints and bounds could not be satisfied.
36         // The problem has no feasible solution.
37         INFORM_NONLINEAR_CONSTRAINTS_INFEASIBLE = 3,
38         // The nonlinear constraints and bounds could not be satisfied.
39         // The problem may have no feasible solution.
40         INFORM_ITERATION_LIMIT = 4,
41         // The major iteration limit was reached (Mx status BLUE).
42         INFORM_NOT_AT_OPTIMUM = 6,
43         // The model does not satisfy the first-order optimality conditions
44         // to the required accuracy, and no improved point for the
45         // merit function could be found during the final linesearch (Mx status RED)
46         INFORM_BAD_DERIVATIVES = 7,
47         // The function derivates returned by funcon or funobj
48         // appear to be incorrect.
49         INFORM_INVALID_PARAM = 9
50         // An input parameter was invalid'
51 };
52
53 struct matrixVectorProdTerm {
54         int hentry;
55         int gentry;
56         int dest;
57         matrixVectorProdTerm() {}
58         matrixVectorProdTerm(int he, int ge, int de) {
59                 hentry=he;
60                 gentry=ge;
61                 dest=de;
62         }
63         bool operator< (const matrixVectorProdTerm &j) const {
64           if (hentry < j.hentry) return true;
65           if (hentry == j.hentry) {
66             if (gentry < j.gentry) return true;
67             if (gentry == j.gentry && dest < j.dest) return true;
68           }
69           return false;
70         };
71         bool operator==(const matrixVectorProdTerm &i) const {
72                 return i.hentry == hentry && i.gentry == gentry && i.dest == dest;
73         }
74 };
75
76 // The idea of FitContext is to eventually enable fitting from
77 // multiple starting values in parallel.
78
79 class FitContext {
80         static omxFitFunction *RFitFunction;
81
82         FitContext *parent;
83
84  public:
85         FreeVarGroup *varGroup;
86         std::vector<int> mapToParent;
87         double mac;
88         double fit;
89         double *est;
90         int *flavor;
91         //      double *denom;
92         double *grad;
93         double *hess;
94         double *ihess;
95         std::vector< matrixVectorProdTerm > hgProd;
96         double caution;
97         int iterations;
98         enum ComputeInform inform;
99         int wanted;
100         bool changedEstimates; // only used for FF_COMPUTE_PREOPTIMIZE
101
102         void init();
103         FitContext(std::vector<double> &startingValues);
104         FitContext(FitContext *parent, FreeVarGroup *group);
105         void fixHessianSymmetry(int want);
106         void copyParamToModel(omxState* os, double *at);
107         void copyParamToModel(omxState *os);
108         void copyParamToModel(omxMatrix *mat, double *at);
109         void copyParamToModel(omxMatrix *mat);
110         double *take(int want);
111         void maybeCopyParamToModel(omxState* os);
112         void updateParent();
113         void updateParentAndFree();
114         void log(const char *where);
115         void log(const char *where, int what);
116         ~FitContext();
117         
118         static void cacheFreeVarDependencies();
119         static void setRFitFunction(omxFitFunction *rff);
120 };
121
122 class omxCompute {
123  public:
124         FreeVarGroup *varGroup;
125         omxCompute();
126         virtual void initFromFrontend(SEXP rObj);
127         virtual omxFitFunction *getFitFunction() { return NULL; }
128         virtual void compute(FitContext *fc) = 0;
129         virtual void reportResults(FitContext *fc, MxRList *out) = 0;
130         virtual double getOptimizerStatus() { return NA_REAL; }  // backward compatibility
131         virtual ~omxCompute();
132 };
133
134 class Ramsay1975 {
135         // Ramsay, J. O. (1975). Solving Implicit Equations in
136         // Psychometric Data Analysis.  Psychometrika, 40(3), 337-360.
137
138         FitContext *fc;
139         size_t numParam;
140         int flavor;
141         int verbose;
142         double minCaution;
143         double highWatermark;
144         std::vector<double> prevEst;
145         std::vector<double> prevAdj1;
146         std::vector<double> prevAdj2;
147
148 public:
149         double maxCaution;
150         double caution;
151
152         Ramsay1975(FitContext *fc, int flavor, double caution, int verbose, double minCaution);
153         void recordEstimate(int px, double newEst);
154         void apply();
155         void recalibrate(bool *restart);
156         void restart();
157 };
158
159 class omxCompute *omxNewCompute(omxState* os, const char *type);
160
161 class omxCompute *newComputeGradientDescent();
162 class omxCompute *newComputeEstimatedHessian();
163 class omxCompute *newComputeNewtonRaphson();
164
165 void omxApproxInvertPosDefTriangular(int dim, double *hess, double *ihess, double *stress);
166 void omxApproxInvertPackedPosDefTriangular(int dim, int *mask, double *packedHess, double *stress);
167
168 #endif