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