Don't hide errors
[openmx:openmx.git] / src / omxNPSOLSpecific.cpp
1 /*
2  *  Copyright 2007-2014 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 #define R_NO_REMAP
19 #include <R.h>
20 #include <Rinternals.h>
21
22 #include "omxState.h"
23 #include "omxNPSOLSpecific.h"
24 #include "omxMatrix.h"
25 #include "glue.h"
26 #include "omxImportFrontendState.h"
27 #include "Compute.h"
28 #include "npsolswitch.h"
29
30 /* NPSOL-specific globals */
31 const double NPSOL_BIGBND = 1e20;
32
33 #if HAS_NPSOL
34 static const char* anonMatrix = "anonymous matrix";
35 static omxMatrix *NPSOL_fitMatrix = NULL;
36 static int NPSOL_currentInterval = -1;
37 static FitContext *NPSOL_fc = NULL;
38 static bool NPSOL_useGradient;
39 static int NPSOL_verbose;
40 #endif
41
42 #ifdef  __cplusplus
43 extern "C" {
44 #endif
45
46 /* NPSOL-related functions */
47 extern void F77_SUB(npsol)(int *n, int *nclin, int *ncnln, int *ldA, int *ldJ, int *ldR, double *A,
48                             double *bl, double *bu, void* funcon, void* funobj, int *inform, int *iter, 
49                             int *istate, double *c, double *cJac, double *clambda, double *f, double *g, double *R,
50                             double *x, int *iw, int *leniw, double *w, int *lenw);
51 extern void F77_SUB(npoptn)(char* string, int Rf_length);
52
53 #ifdef  __cplusplus
54 }
55 #endif
56
57 // NOTE: All non-linear constraints are applied regardless of free TODO
58 // variable group.  This is probably wrong. TODO
59 void omxSetupBoundsAndConstraints(FreeVarGroup *freeVarGroup, double * bl, double * bu)
60 {
61         size_t n = freeVarGroup->vars.size();
62
63         /* Set min and max limits */
64         for(size_t index = 0; index < n; index++) {
65                 bl[index] = freeVarGroup->vars[index]->lbound;
66                 bu[index] = freeVarGroup->vars[index]->ubound;
67         }
68
69         int index = n;
70         for(int constraintIndex = 0; constraintIndex < globalState->numConstraints; constraintIndex++) {                // Nonlinear constraints:
71                 if(OMX_DEBUG) { mxLog("Constraint %d: ", constraintIndex);}
72                 switch(globalState->conList[constraintIndex].opCode) {
73                 case 0:                                                                 // Less than: Must be strictly less than 0.
74                         if(OMX_DEBUG) { mxLog("Bounded at (-INF, 0).");}
75                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
76                                 bl[index] = NEG_INF;
77                                 bu[index] = -0.0;
78                                 index++;
79                         }
80                         break;
81                 case 1:                                                                 // Equal: Must be roughly equal to 0.
82                         if(OMX_DEBUG) { mxLog("Bounded at (-0, 0).");}
83                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
84                                 bl[index] = -0.0;
85                                 bu[index] = 0.0;
86                                 index++;
87                         }
88                         break;
89                 case 2:                                                                 // Greater than: Must be strictly greater than 0.
90                         if(OMX_DEBUG) { mxLog("Bounded at (0, INF).");}
91                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
92                                 if(OMX_DEBUG) { mxLog("\tBounds set for constraint %d.%d.", constraintIndex, offset);}
93                                 bl[index] = 0.0;
94                                 bu[index] = INF;
95                                 index++;
96                         }
97                         break;
98                 default:
99                         if(OMX_DEBUG) { mxLog("Bounded at (-INF, INF).");}
100                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
101                                 bl[index] = NEG_INF;
102                                 bu[index] = INF;
103                                 index++;
104                         }
105                         break;
106                 }
107         }
108 }
109
110 #if HAS_NPSOL
111
112 /****** Objective Function *********/
113 static void
114 npsolObjectiveFunction1(int* mode, int* n, double* x,
115                         double* f, double* g, double *hessian, int* nstate )
116 {
117         unsigned short int checkpointNow = FALSE;
118
119         if(OMX_DEBUG) {mxLog("Starting Objective Run.");}
120
121         if(*mode == 1) {
122                 omxSetMajorIteration(globalState, globalState->majorIteration + 1);
123                 omxSetMinorIteration(globalState, 0);
124                 checkpointNow = TRUE;                                   // Only checkpoint at major iterations.
125         } else omxSetMinorIteration(globalState, globalState->minorIteration + 1);
126
127         omxMatrix* fitMatrix = NPSOL_fitMatrix;
128
129         /* Interruptible? */
130         R_CheckUserInterrupt();
131
132         NPSOL_fc->copyParamToModel(globalState, x);
133
134         if (*mode > 0 && NPSOL_useGradient &&
135             fitMatrix->fitFunction->gradientAvailable && NPSOL_currentInterval < 0) {
136                 size_t numParams = NPSOL_fc->varGroup->vars.size();
137                 NPSOL_fc->grad = Eigen::VectorXd::Zero(numParams);
138
139                 omxFitFunctionCompute(fitMatrix->fitFunction, FF_COMPUTE_FIT|FF_COMPUTE_GRADIENT, NPSOL_fc);
140                 if (NPSOL_verbose) {
141                         NPSOL_fc->log(FF_COMPUTE_FIT|FF_COMPUTE_ESTIMATE|FF_COMPUTE_GRADIENT);
142                 }
143         } else {
144                 omxFitFunctionCompute(fitMatrix->fitFunction, FF_COMPUTE_FIT, NPSOL_fc);
145                 if (NPSOL_verbose) {
146                         NPSOL_fc->log(FF_COMPUTE_FIT|FF_COMPUTE_ESTIMATE);
147                 }
148         }
149
150         if (!R_FINITE(fitMatrix->data[0])) {
151                 *mode = -1;
152         } else {
153                 NPSOL_fc->resetIterationError();
154         }
155
156         *f = fitMatrix->data[0];
157         if(OMX_VERBOSE) {
158                 mxLog("Fit function value is: %f, Mode is %d.", fitMatrix->data[0], *mode);
159         }
160
161         if(OMX_DEBUG) { mxLog("-======================================================-"); }
162
163         if(checkpointNow && globalState->numCheckpoints != 0) { // If it's a new major iteration
164                 omxSaveCheckpoint(x, *f, FALSE);                // Check about saving a checkpoint
165         }
166
167 }
168
169 void F77_SUB(npsolObjectiveFunction)
170         (       int* mode, int* n, double* x,
171                 double* f, double* g, int* nstate )
172 {
173         npsolObjectiveFunction1(mode, n, x, f, g, NULL, nstate);
174 }
175
176 /* Objective function for confidence interval limit finding. 
177  * Replaces the standard objective function when finding confidence intervals. */
178 void F77_SUB(npsolLimitObjectiveFunction)
179         (       int* mode, int* n, double* x, double* f, double* g, int* nstate ) {
180                 
181                 if(OMX_VERBOSE) mxLog("Calculating interval %d, %s boundary:", NPSOL_currentInterval, (Global->intervalList[NPSOL_currentInterval].calcLower?"lower":"upper"));
182
183                 F77_CALL(npsolObjectiveFunction)(mode, n, x, f, g, nstate);     // Standard objective function call
184
185                 omxConfidenceInterval *oCI = &(Global->intervalList[NPSOL_currentInterval]);
186                 
187                 omxRecompute(oCI->matrix);
188                 
189                 double CIElement = omxMatrixElement(oCI->matrix, oCI->row, oCI->col);
190
191                 if(OMX_DEBUG) {
192                         mxLog("Finding Confidence Interval Likelihoods: lbound is %f, ubound is %f, estimate likelihood is %f, and element current value is %f.",
193                                 oCI->lbound, oCI->ubound, *f, CIElement);
194                 }
195
196                 /* Catch boundary-passing condition */
197                 if(std::isnan(CIElement) || std::isinf(CIElement)) {
198                         NPSOL_fc->recordIterationError("Confidence interval is in a range that is currently incalculable. Add constraints to keep the value in the region where it can be calculated.");
199                         *mode = -1;
200                         return;
201                 }
202
203                 if(oCI->calcLower) {
204                         double diff = oCI->lbound - *f;         // Offset - likelihood
205                         *f = diff * diff + CIElement;
206                                 // Minimize element for lower bound.
207                 } else {
208                         double diff = oCI->ubound - *f;                 // Offset - likelihood
209                         *f = diff * diff - CIElement;
210                                 // Maximize element for upper bound.
211                 }
212
213                 if(OMX_DEBUG) {
214                         mxLog("Interval fit function in previous iteration was calculated to be %f.", *f);
215                 }
216 }
217
218 /* (Non)Linear Constraint Functions */
219 void F77_SUB(npsolConstraintFunction)
220         (       int *mode, int *ncnln, int *n,
221                 int *ldJ, int *needc, double *x,
222                 double *c, double *cJac, int *nstate)
223 {
224
225         if(OMX_DEBUG) { mxLog("Constraint function called.");}
226
227         if(*mode==1) {
228                 if(OMX_DEBUG) {
229                         mxLog("But only gradients requested.  Returning.");
230                         mxLog("-=====================================================-");
231                 }
232
233                 return;
234         }
235
236         int j, k, l = 0;
237
238         NPSOL_fc->copyParamToModel(globalState, x);
239
240         for(j = 0; j < globalState->numConstraints; j++) {
241                 omxRecompute(globalState->conList[j].result);
242                 if(OMX_VERBOSE) { omxPrint(globalState->conList[j].result, "Constraint evaluates as:"); }
243                 for(k = 0; k < globalState->conList[j].size; k++){
244                         c[l++] = globalState->conList[j].result->data[k];
245                 }
246         }
247         if(OMX_DEBUG) { mxLog("-=======================================================-"); }
248 }
249
250 void omxInvokeNPSOL(omxMatrix *fitMatrix, FitContext *fc,
251                     int *inform_out, int *iter_out, bool useGradient, FreeVarGroup *freeVarGroup,
252                     int verbose, double *hessOut)
253 {
254         // Will fail if we re-enter after an exception
255         //if (NPSOL_fitMatrix) Rf_error("NPSOL is not reentrant");
256         NPSOL_fitMatrix = fitMatrix;
257         NPSOL_verbose = verbose;
258
259         NPSOL_useGradient = useGradient;
260         NPSOL_fc = fc;
261         double *x = fc->est;
262         fc->grad.resize(fc->numParam); // ensure memory is allocated
263         double *g = fc->grad.data();
264
265     double *A=NULL, *bl=NULL, *bu=NULL, *c=NULL, *clambda=NULL, *w=NULL; //  *g, *R, *cJac,
266  
267     int k, ldA, ldJ, ldR, inform, iter, leniw, lenw; 
268  
269     double *cJac = NULL;    // Hessian (Approx) and Jacobian
270  
271     int *iw = NULL;
272  
273     int *istate = NULL;                 // Current state of constraints (0 = no, 1 = lower, 2 = upper, 3 = both (equality))
274  
275     int nctotl, nlinwid, nlnwid;    // Helpful side variables.
276  
277     int nclin = 0;
278     int ncnln = globalState->ncnln;
279  
280     /* NPSOL Arguments */
281     void (*funcon)(int*, int*, int*, int*, int*, double*, double*, double*, int*);
282  
283     funcon = F77_SUB(npsolConstraintFunction);
284  
285         /* Set boundaries and widths. */
286         if(nclin <= 0) {
287             nclin = 0;                  // This should never matter--nclin should always be non-negative.
288             nlinwid = 1;                // For memory allocation purposes, nlinwid > 0
289         } else {                        // nlinwid is  used to calculate ldA, and eventually the size of A.
290             nlinwid = nclin;
291         }
292  
293         if(ncnln <= 0) {
294             ncnln = 0;                  // This should never matter--ncnln should always be non-negative.
295             nlnwid = 1;                 // For memory allocation purposes nlnwid > 0
296         } else {                        // nlnwid is used to calculate ldJ, and eventually the size of J.
297             nlnwid = ncnln;
298         }
299  
300         int n = int(freeVarGroup->vars.size());
301  
302         nctotl = n + nlinwid + nlnwid;
303  
304         leniw = 3 * n + nclin + 2 * ncnln;
305         lenw = 2 * n * n + n * nclin + 2 * n * ncnln + 20 * n + 11 * nclin + 21 * ncnln;
306  
307         ldA = nlinwid;          // NPSOL specifies this should be either 1 or nclin, whichever is greater
308         ldJ = nlnwid;           // NPSOL specifies this should be either 1 or nclin, whichever is greater
309         ldR = n;                // TODO: Test alternative versions of the size of R to see what's best.
310  
311     /* Allocate arrays */
312         A       = (double*) R_alloc (ldA * n, sizeof ( double )  );
313         bl      = (double*) R_alloc ( nctotl, sizeof ( double ) );
314         bu      = (double*) R_alloc (nctotl, sizeof ( double ) );
315         c       = (double*) R_alloc (nlnwid, sizeof ( double ));
316         cJac    = (double*) R_alloc (ldJ * n, sizeof ( double ) );
317         clambda = (double*) R_alloc (nctotl, sizeof ( double )  );
318         w       = (double*) R_alloc (lenw, sizeof ( double ));
319         istate  = (int*) R_alloc (nctotl, sizeof ( int ) );
320         iw      = (int*) R_alloc (leniw, sizeof ( int ));
321  
322         /* Set up actual run */
323  
324         omxSetupBoundsAndConstraints(freeVarGroup, bl, bu);
325  
326         /* Initialize Starting Values */
327         if(OMX_VERBOSE) {
328             mxLog("--------------------------");
329             mxLog("Starting Values (%d) are:", n);
330         }
331         for(k = 0; k < n; k++) {
332             if((x[k] == 0.0)) {
333                 x[k] += 0.1;
334             }
335             if(OMX_VERBOSE) { mxLog("%d: %f", k, x[k]); }
336         }
337         if(OMX_DEBUG) {
338             mxLog("--------------------------");
339             mxLog("Setting up optimizer...");
340         }
341  
342     /*  F77_CALL(npsol)
343         (   int *n,                 -- Number of variables
344             int *nclin,             -- Number of linear constraints
345             int *ncnln,             -- Number of nonlinear constraints
346             int *ldA,               -- Row dimension of A (Linear Constraints)
347             int *ldJ,               -- Row dimension of cJac (Jacobian)
348             int *ldR,               -- Row dimension of R (Hessian)
349             double *A,              -- Linear Constraints Array A (in Column-major order)
350             double *bl,             -- Lower Bounds Array (at least n + nclin + ncnln long)
351             double *bu,             -- Upper Bounds Array (at least n + nclin + ncnln long)
352             function funcon,        -- Nonlinear constraint function
353             function funobj,        -- Objective function
354             int *inform,            -- Used to report state.  Need not be initialized.
355             int *iter,              -- Used to report number of major iterations performed.  Need not be initialized.
356             int *istate,            -- Initial State.  Need not be initialized unless using Warm Start.
357             double *c,              -- Array of Rf_length ncnln.  Need not be initialized.  Reports nonlinear constraints at final iteration.
358             double *cJac,           -- Array of Row-Rf_length ldJ.  Unused if ncnln = 0. Generally need not be initialized.
359             double *clambda,        -- Array of Rf_length n+nclin+ncnln.  Need not be initialized unless using Warm Start. Reports final QP multipliers.
360             double *f,              -- Used to report final objective value.  Need not be initialized.
361             double *g,              -- Array of Rf_length n. Used to report final objective gradient.  Need not be initialized.
362             double *R,              -- Array of Rf_length ldR.  Need not be intialized unless using Warm Start.
363             double *x,              -- Array of Rf_length n.  Contains initial solution estimate.
364             int *iw,                -- Array of Rf_length leniw. Need not be initialized.  Provides workspace.
365             int *leniw,             -- Length of iw.  Must be at least 3n + nclin + ncnln.
366             double *w,              -- Array of Rf_length lenw. Need not be initialized.  Provides workspace.
367             int *lenw               -- Length of w.  Must be at least 2n^2 + n*nclin + 2*n*ncnln + 20*n + 11*nclin +21*ncnln
368         )
369  
370         bl, bu, istate, and clambda are all Rf_length n+nclin+ncnln.
371             First n elements refer to the vars, in order.
372             Next nclin elements refer to bounds on Ax
373             Last ncnln elements refer to bounds on c(x)
374  
375         All arrays must be in column-major order.
376  
377         */
378  
379         if(OMX_DEBUG) {
380             mxLog("Set.");
381         }
382  
383         F77_CALL(npsol)(&n, &nclin, &ncnln, &ldA, &ldJ, &ldR, A, bl, bu, (void*)funcon,
384                         (void*) F77_SUB(npsolObjectiveFunction), &inform, &iter, istate, c, cJac,
385                         clambda, &fc->fit, g, hessOut, x, iw, &leniw, w, &lenw);
386
387         if(OMX_DEBUG) { mxLog("Final Objective Value is: %f", fc->fit); }
388  
389         omxSaveCheckpoint(x, fc->fit, TRUE);
390  
391         NPSOL_fc->copyParamToModel(globalState);
392  
393     *inform_out = inform;
394     *iter_out   = iter;
395  
396     NPSOL_fitMatrix = NULL;
397     NPSOL_fc = NULL;
398 }
399  
400  
401 // Mostly duplicated code in omxCSOLNPConfidenceIntervals
402 // needs to be refactored so there is only 1 copy of CI
403 // code that can use whatever optimizer is provided.
404 void omxNPSOLConfidenceIntervals(omxMatrix *fitMatrix, FitContext *fc)
405 {
406         int ciMaxIterations = Global->ciMaxIterations;
407         // Will fail if we re-enter after an exception
408         //if (NPSOL_fitMatrix) Rf_error("NPSOL is not reentrant");
409         NPSOL_fitMatrix = fitMatrix;
410         NPSOL_fc = fc;
411         FreeVarGroup *freeVarGroup = fitMatrix->fitFunction->freeVarGroup;
412
413     double *A=NULL, *bl=NULL, *bu=NULL, *c=NULL, *clambda=NULL, *w=NULL; //  *g, *R, *cJac,
414  
415     int ldA, ldJ, ldR, inform, iter, leniw, lenw; 
416  
417     double *cJac = NULL;    // Hessian (Approx) and Jacobian
418  
419     int *iw = NULL;
420  
421     int *istate = NULL;                 // Current state of constraints (0 = no, 1 = lower, 2 = upper, 3 = both (equality))
422  
423     int nctotl, nlinwid, nlnwid;    // Helpful side variables.
424  
425     int n = int(freeVarGroup->vars.size());
426     double optimum = fc->fit;
427     double *optimalValues = fc->est;
428     double f = optimum;
429     std::vector< double > x(n, *optimalValues);
430     std::vector< double > gradient(n);
431     std::vector< double > hessian(n * n);
432
433     /* NPSOL Arguments */
434     void (*funcon)(int*, int*, int*, int*, int*, double*, double*, double*, int*);
435  
436     funcon = F77_SUB(npsolConstraintFunction);
437  
438     int nclin = 0;
439     int ncnln = globalState->ncnln;
440  
441     /* Set boundaries and widths. */
442     if(nclin <= 0) {
443         nclin = 0;                  // This should never matter--nclin should always be non-negative.
444         nlinwid = 1;                // For memory allocation purposes, nlinwid > 0
445     } else {                        // nlinwid is  used to calculate ldA, and eventually the size of A.
446         nlinwid = nclin;
447     }
448  
449     if(ncnln <= 0) {
450         ncnln = 0;                  // This should never matter--ncnln should always be non-negative.
451         nlnwid = 1;                 // For memory allocation purposes nlnwid > 0
452     } else {                        // nlnwid is used to calculate ldJ, and eventually the size of J.
453         nlnwid = ncnln;
454     }
455  
456     nctotl = n + nlinwid + nlnwid;
457  
458     leniw = 3 * n + nclin + 2 * ncnln;
459     lenw = 2 * n * n + n * nclin + 2 * n * ncnln + 20 * n + 11 * nclin + 21 * ncnln;
460  
461     ldA = nlinwid;          // NPSOL specifies this should be either 1 or nclin, whichever is greater
462     ldJ = nlnwid;           // NPSOL specifies this should be either 1 or nclin, whichever is greater
463     ldR = n;                // TODO: Test alternative versions of the size of R to see what's best.
464  
465     /* Allocate arrays */
466     A       = (double*) R_alloc (ldA * n, sizeof ( double )  );
467     bl      = (double*) R_alloc ( nctotl, sizeof ( double ) );
468     bu      = (double*) R_alloc (nctotl, sizeof ( double ) );
469     c       = (double*) R_alloc (nlnwid, sizeof ( double ));
470     cJac    = (double*) R_alloc (ldJ * n, sizeof ( double ) );
471     clambda = (double*) R_alloc (nctotl, sizeof ( double )  );
472     w       = (double*) R_alloc (lenw, sizeof ( double ));
473     istate  = (int*) R_alloc (nctotl, sizeof ( int ) );
474     iw      = (int*) R_alloc (leniw, sizeof ( int ));
475  
476  
477     omxSetupBoundsAndConstraints(freeVarGroup, bl, bu);
478  
479         if(OMX_DEBUG) { mxLog("Calculating likelihood-based confidence intervals."); }
480
481         for(int i = 0; i < Global->numIntervals; i++) {
482
483                         omxConfidenceInterval *currentCI = &(Global->intervalList[i]);
484
485                         int msgLength = 45;
486  
487                         if (currentCI->matrix->name == NULL) {
488                                 msgLength += strlen(anonMatrix);
489                         } else {
490                                 msgLength += strlen(currentCI->matrix->name);
491                         }
492             
493                         char *message = Calloc(msgLength, char);
494  
495                         if (currentCI->matrix->name == NULL) {
496                                 snprintf(message, msgLength, "%s[%d, %d] begin lower interval",
497                                         anonMatrix, currentCI->row + 1, currentCI->col + 1);
498                         } else {
499                                 snprintf(message, msgLength, "%s[%d, %d] begin lower interval",
500                                         currentCI->matrix->name, currentCI->row + 1, currentCI->col + 1);
501                         }
502  
503                         omxWriteCheckpointMessage(message);
504  
505                         memcpy(x.data(), optimalValues, n * sizeof(double)); // Reset to previous optimum
506                         NPSOL_currentInterval = i;
507
508             currentCI->lbound += optimum;          // Convert from offsets to targets
509             currentCI->ubound += optimum;          // Convert from offsets to targets
510  
511             /* Set up for the lower bound */
512             inform = -1;
513             // Number of times to keep trying.
514             int cycles = ciMaxIterations;
515             double value = INF;
516             double objDiff = 1.e-4;     // TODO : Use function precision to determine CI jitter?
517             while(inform != 0 && cycles > 0) {
518                 /* Find lower limit */
519                 currentCI->calcLower = TRUE;
520                 F77_CALL(npsol)(&n, &nclin, &ncnln, &ldA, &ldJ, &ldR, A, bl, bu, (void*)funcon,
521                     (void*) F77_SUB(npsolLimitObjectiveFunction), &inform, &iter, istate, c, cJac,
522                                 clambda, &f, gradient.data(), hessian.data(), x.data(), iw, &leniw, w, &lenw);
523  
524                 currentCI->lCode = inform;
525                 if(f < value) {
526                     currentCI->min = omxMatrixElement(currentCI->matrix, currentCI->row, currentCI->col);
527                     value = f;
528                     omxSaveCheckpoint(x.data(), f, TRUE);
529                 }
530  
531                 if(inform != 0 && OMX_DEBUG) {
532                     mxLog("Calculation of lower interval %d failed: Bad inform value of %d",
533                             i, inform);
534                 }
535                 cycles--;
536                 if(inform != 0) {
537                     unsigned int jitter = TRUE;
538                     for(int j = 0; j < n; j++) {
539                         if(fabs(x[j] - optimalValues[j]) > objDiff) {
540                             jitter = FALSE;
541                             break;
542                         }
543                     }
544                     if(jitter) {
545                         for(int j = 0; j < n; j++) {
546                             x[j] = optimalValues[j] + objDiff;
547                         }
548                     }
549                 }
550             }
551  
552             if(OMX_DEBUG) { mxLog("Found lower bound %d.  Seeking upper.", i); }
553             // TODO: Repopulate original optimizer state in between CI calculations
554
555                         if (currentCI->matrix->name == NULL) {
556                                 snprintf(message, msgLength, "%s[%d, %d] begin upper interval", 
557                                         anonMatrix, currentCI->row + 1, currentCI->col + 1);
558                         } else {
559                                 snprintf(message, msgLength, "%s[%d, %d] begin upper interval",
560                                         currentCI->matrix->name, currentCI->row + 1, currentCI->col + 1);
561                         }
562  
563                         omxWriteCheckpointMessage(message);
564  
565                         Free(message);
566  
567                         memcpy(x.data(), optimalValues, n * sizeof(double));
568  
569             /* Reset for the upper bound */
570             value = INF;
571             inform = -1;
572             cycles = ciMaxIterations;
573  
574             while(inform != 0 && cycles >= 0) {
575                 /* Find upper limit */
576                 currentCI->calcLower = FALSE;
577                 F77_CALL(npsol)(&n, &nclin, &ncnln, &ldA, &ldJ, &ldR, A, bl, bu, (void*)funcon,
578                                     (void*) F77_SUB(npsolLimitObjectiveFunction), &inform, &iter, istate, c, cJac,
579                                 clambda, &f, gradient.data(), hessian.data(), x.data(), iw, &leniw, w, &lenw);
580  
581                 currentCI->uCode = inform;
582                 if(f < value) {
583                     currentCI->max = omxMatrixElement(currentCI->matrix, currentCI->row, currentCI->col);
584                     value = f;
585                     omxSaveCheckpoint(x.data(), f, TRUE);
586                 }
587  
588                 if(inform != 0 && OMX_DEBUG) {
589                     mxLog("Calculation of upper interval %d failed: Bad inform value of %d",
590                             i, inform);
591                 }
592                 cycles--;
593                 if(inform != 0) {
594                     unsigned int jitter = TRUE;
595                     for(int j = 0; j < n; j++) {
596                         if(fabs(x[j] - optimalValues[j]) > objDiff){
597                             jitter = FALSE;
598                             break;
599                         }
600                     }
601                     if(jitter) {
602                         for(int j = 0; j < n; j++) {
603                             x[j] = optimalValues[j] + objDiff;
604                         }
605                     }
606                 }
607             }
608             if(OMX_DEBUG) {mxLog("Found Upper bound %d.", i);}
609         }
610
611         NPSOL_fc = NULL;
612         NPSOL_fitMatrix = NULL;
613         NPSOL_currentInterval = -1;
614 }
615  
616 void omxSetNPSOLOpts(SEXP options)
617 {
618         static const char *whitelist[] = {
619                 "Central Difference Interval",
620                 "Crash Tolerance",
621                 "Derivative level",
622                 "Difference interval",
623                 "Feasibility tolerance",
624                 "Function precision",
625                 "Hessian",
626                 "Infinite bound size",
627                 "Infinite step size",
628                 "Iteration limit",
629                 "Iters",
630                 "Line search tolerance",
631                 "Major iteration limit",
632                 "Major iterations",
633                 "Print level",
634                 "Print file",
635                 "Minor iteration limit",
636                 "Minor print level",
637                 "Nolist",
638                 "Optimality tolerance",
639                 "Step limit",
640                 "Summary file",
641                 "Verify level",
642                 0
643         };
644
645         const int opBufLen = 250;
646         char optionCharArray[opBufLen];
647         int numOptions = Rf_length(options);
648         SEXP optionNames;
649         Rf_protect(optionNames = Rf_getAttrib(options, R_NamesSymbol));
650                 for(int i = 0; i < numOptions; i++) {
651                         const char *nextOptionName = CHAR(STRING_ELT(optionNames, i));
652                         const char *nextOptionValue = CHAR(Rf_asChar(VECTOR_ELT(options, i)));
653                         bool ok=false;
654                         for (int wx=0; whitelist[wx]; ++wx) {
655                                 if (matchCaseInsensitive(nextOptionName, whitelist[wx])) {
656                                         ok=true;
657                                         break;
658                                 }
659                         }
660                         if (!ok) continue;
661                         snprintf(optionCharArray, opBufLen, "%s %s", nextOptionName, nextOptionValue);
662                         F77_CALL(npoptn)(optionCharArray, strlen(optionCharArray));
663                         if(OMX_DEBUG) { mxLog("Option %s ", optionCharArray); }
664                 }
665                 Rf_unprotect(1); // optionNames
666 }
667
668 #endif