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