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