Csolnp_warnings removed v2
[openmx:openmx.git] / src / omxCsolnp.cpp
1 /*
2  *  Copyright 2007-2012 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 #include "omxState.h"
22 #include "omxNPSOLSpecific.h"
23 #include "omxMatrix.h"
24 #include "glue.h"
25 #include "omxImportFrontendState.h"
26 #include "matrix.h"
27 #include "omxCsolnp.h"
28
29 static const char* anonMatrix = "anonymous matrix";
30
31 /* NPSOL-related functions */
32 //************************* npsol ****************************//
33 //int solnp(Matrix solPars, double (*solFun)(Matrix),  Matrix solEqB,  Matrix (*solEqBFun)( Matrix),  Matrix (*solEqBStartFun)(Matrix),  Matrix solLB,  Matrix solUB,  Matrix solIneqUB,  Matrix solIneqLB,  Matrix solctrl, bool debugToggle);
34
35 static omxMatrix *GLOB_fitMatrix = NULL;
36 static FitContext *GLOB_fc = NULL;
37 static int CSOLNP_currentInterval = -1;
38
39 Matrix fillMatrix(int cols, int rows, double* array)
40 {
41     Matrix t = new_matrix(cols, rows);
42         int i,j;
43         for(i=0;i<rows;i++){
44                 for(j=0;j<cols;j++) {
45                         M(t,j,i)=array[j];
46                 }
47         }
48         return t;
49 }
50
51
52 //****** Objective Function *********//
53 double csolnpObjectiveFunction(Matrix myPars, int verbose)
54 {
55             
56         unsigned short int checkpointNow = FALSE;
57     
58         if(OMX_DEBUG) {printf("Starting Objective Run.");}
59     
60         omxMatrix* fitMatrix = GLOB_fitMatrix;
61     
62         omxResetStatus(globalState);                                            // Clear Error State recursively
63     
64     /* Interruptible? */
65         R_CheckUserInterrupt();
66     /* This allows for abitrary repopulation of the free parameters.
67      * Typically, the default is for repopulateFun to be NULL,
68      * and then handleFreeVarList is invoked */
69     
70         GLOB_fc->copyParamToModel(globalState, myPars.t);
71     
72     
73     omxFitFunctionCompute(fitMatrix->fitFunction, FF_COMPUTE_FIT, GLOB_fc);
74
75     
76     if (!R_FINITE(fitMatrix->data[0])) {
77         omxRaiseErrorf(globalState, "Fit function returned %g at iteration %d.%d",
78                        fitMatrix->data[0], globalState->majorIteration, globalState->minorIteration);
79     }
80     if(isErrorRaised(globalState))
81     {
82         if(OMX_DEBUG)
83         {    mxLog("Error status reported.");}
84     }
85     if (isErrorRaised(globalState) || !R_FINITE(fitMatrix->data[0]))
86     {    GLOB_fc->fit = 1e24;}
87     else {  GLOB_fc->fit = fitMatrix->data[0];}
88     
89     
90         if(verbose >= 1) {
91                 mxLog("Fit function value is: %.32f\n", fitMatrix->data[0]);
92         }
93     
94         if(checkpointNow && globalState->numCheckpoints != 0) { // If it's a new major iteration
95                 omxSaveCheckpoint(myPars.t, GLOB_fc->fit, FALSE);               // Check about saving a checkpoint
96         }
97         return GLOB_fc->fit;
98 }
99
100
101 /* Objective function for confidence interval limit finding.
102  * Replaces the standard objective function when finding confidence intervals. */
103 double csolnpLimitObjectiveFunction(Matrix myPars, int verbose)
104 {
105     //double* f = NULL;
106         if (verbose >= 3) {
107                 printf("myPars inside obj is: ");
108                 print(myPars); putchar('\n');
109         }
110     
111     if(OMX_VERBOSE) mxLog("Calculating interval %d, %s boundary:", CSOLNP_currentInterval, (Global->intervalList[CSOLNP_currentInterval].calcLower?"lower":"upper"));
112     
113     GLOB_fc->fit = csolnpObjectiveFunction(myPars, verbose);
114     
115     omxConfidenceInterval *oCI = &(Global->intervalList[CSOLNP_currentInterval]);
116     
117     omxRecompute(oCI->matrix);
118     
119     double CIElement = omxMatrixElement(oCI->matrix, oCI->row, oCI->col);
120     
121     if(OMX_DEBUG) {
122         mxLog("Finding Confidence Interval Likelihoods: lbound is %f, ubound is %f, estimate likelihood is %f, and element current value is %f.",
123               oCI->lbound, oCI->ubound, GLOB_fc->fit, CIElement);
124     }
125     
126     /* Catch boundary-passing condition */
127     if(isnan(CIElement) || isinf(CIElement)) {
128         omxRaiseError(globalState, -1,
129                       "Confidence interval is in a range that is currently incalculable. Add constraints to keep the value in the region where it can be calculated.");
130         return GLOB_fc->fit;
131     }
132     
133     if(oCI->calcLower) {
134         double diff = oCI->lbound - GLOB_fc->fit;               // Offset - likelihood
135         GLOB_fc->fit = diff * diff + CIElement;
136         // Minimize element for lower bound.
137     } else {
138         double diff = oCI->ubound - GLOB_fc->fit;                       // Offset - likelihood
139         GLOB_fc->fit = diff * diff - CIElement;
140         // Maximize element for upper bound.
141     }
142     
143     if(OMX_DEBUG) {
144         mxLog("Interval fit function in previous iteration was calculated to be %f.", GLOB_fc->fit);
145     }
146     
147     return GLOB_fc->fit;
148 }
149
150
151 /* (Non)Linear Constraint Functions */
152 Matrix csolnpEqualityFunction(int verbose)
153 {
154         int j, k, eq_n = 0;
155     int l = 0;
156     double EMPTY = -999999.0;
157     Matrix myEqBFun;
158     
159     if (verbose) mxLog("Starting csolnpEqualityFunction.");
160     
161     for(j = 0; j < globalState->numConstraints; j++) {
162             if (globalState->conList[j].opCode == 1) {
163                     eq_n += globalState->conList[j].size;
164             }
165     }
166     
167     if (verbose >= 1) {
168             mxLog("no.of constraints is: %d.", globalState->numConstraints);
169             mxLog("neq is: %d.", eq_n);
170     }
171     
172     if (eq_n == 0)
173     {
174         myEqBFun = fill(1, 1, EMPTY);
175     }
176     else {
177             myEqBFun = fill(eq_n, 1, EMPTY);
178             
179             for(j = 0; j < globalState->numConstraints; j++) {
180                     if (globalState->conList[j].opCode == 1) {
181                             if (verbose >= 1) {
182                                     mxLog("result is: %2f", globalState->conList[j].result->data[0]);
183                             }
184                             omxRecompute(globalState->conList[j].result);
185                             if (verbose >= 1) {
186                                     mxLog("%.16f", globalState->conList[j].result->data[0]);
187                                     mxLog("size is: %d", globalState->conList[j].size);
188                             }
189                     }
190                     for(k = 0; k < globalState->conList[j].size; k++){
191                             M(myEqBFun,l,0) = globalState->conList[j].result->data[k];
192                             l = l + 1;
193                     }
194             }
195     }
196     if (verbose >= 1) {
197             printf("myEqBFun is: ");
198             print(myEqBFun);
199             putchar('\n');
200     }
201     return myEqBFun;
202 }
203
204
205 Matrix csolnpIneqFun(int verbose)
206 {
207         int j, k, ineq_n = 0;
208     int l = 0;
209     double EMPTY = -999999.0;
210     Matrix myIneqFun;
211     
212     if (verbose) mxLog("Starting csolnpIneqFun.");
213         
214         for(j = 0; j < globalState->numConstraints; j++) {
215                 if ((globalState->conList[j].opCode == 0) || (globalState->conList[j].opCode == 2))
216         {
217             ineq_n += globalState->conList[j].size;
218         }
219     }
220     
221         if (verbose >= 1) {
222                 printf("no.of constraints is: %d.", globalState->numConstraints); putchar('\n');
223                 printf("ineq_n is: %d.", ineq_n); putchar('\n');
224         }
225     
226     if (ineq_n == 0)
227     {
228         myIneqFun = fill(1, 1, EMPTY);
229     }
230     else
231     {
232         myIneqFun = fill(ineq_n, 1, EMPTY);
233         
234         for(j = 0; j < globalState->numConstraints; j++) {
235             if ((globalState->conList[j].opCode == 0) || globalState->conList[j].opCode == 2)
236             {   omxRecompute(globalState->conList[j].result);}
237             for(k = 0; k < globalState->conList[j].size; k++){
238                 M(myIneqFun,l,0) = globalState->conList[j].result->data[k];
239                 l = l + 1;
240                 
241             }
242         }
243     }
244     
245     return myIneqFun;
246 }
247
248 void omxInvokeCSOLNP(omxMatrix *fitMatrix, FitContext *fc,
249                      int *inform_out, int *iter_out, FreeVarGroup *freeVarGroup,
250                      int verbose)
251
252 {
253         freeMatrices(); // maybe left overs from an aborted optimization attempt
254     
255         GLOB_fitMatrix = fitMatrix;
256         GLOB_fc = fc;
257     
258     double *x = fc->est;
259     //double *g = fc->grad;
260     
261     
262     int k, iter = -1;
263     int inform = 0;
264     
265     double *bl=NULL, *bu=NULL;
266     
267     //double *cJac = NULL;    // Hessian (Approx) and Jacobian
268     
269     int ncnln = globalState->ncnln;
270     int n = int(freeVarGroup->vars.size());
271     
272     double EMPTY = -999999.0;
273     
274     Param_Obj p_obj;
275     Matrix param_hess;
276     Matrix myhess = fill(n*n, 1, (double)0.0);
277     Matrix mygrad;
278     Matrix solIneqLB;
279     Matrix solIneqUB;
280     Matrix solEqB;
281     
282     Matrix myPars = fillMatrix(n, 1, fc->est);
283     
284     double (*solFun)(struct Matrix myPars, int verbose);
285     solFun = &csolnpObjectiveFunction;
286     Matrix (*solEqBFun)(int verbose);
287     solEqBFun = &csolnpEqualityFunction;
288     Matrix (*solIneqFun)(int verbose);
289     solIneqFun = &csolnpIneqFun;
290     
291     /* Set boundaries and widths. */
292     
293     /* Allocate arrays */
294     bl      = Calloc ( n + ncnln, double );
295     bu      = Calloc ( n + ncnln, double );
296     
297     struct Matrix myControl = fill(6,1,(double)0.0);
298     M(myControl,0,0) = 1.0;
299     M(myControl,1,0) = 400.0;
300     M(myControl,2,0) = 800.0;
301     M(myControl,3,0) = 1.0e-7;
302     M(myControl,4,0) = 1.0e-8;
303     M(myControl,5,0) = 0.0;
304     
305     bool myDEBUG = false;
306     /* Set up actual run */
307     
308     /* needs treatment*/
309     if (ncnln == 0)
310     {
311         solIneqLB = fill(1, 1, EMPTY);
312         solIneqUB = fill(1, 1, EMPTY);
313         solEqB = fill(1, 1, EMPTY);
314     }
315     else{
316         int j;
317         int nineqn;
318         int eqn = 0;
319         for(j = 0; j < globalState->numConstraints; j++) {
320             if (globalState->conList[j].opCode == 1)
321             {
322                 eqn += globalState->conList[j].size;
323             }
324         }
325         if (eqn == ncnln) nineqn = 1;
326         else nineqn = ncnln - eqn;
327         
328         solIneqLB = fill(nineqn, 1, EMPTY);
329         solIneqUB = fill(nineqn, 1, EMPTY);
330             if (eqn == 0) {
331                     solEqB = fill(1, 1, EMPTY);
332             } else {
333                     solEqB = fill(eqn, 1, EMPTY);
334             }
335         
336         omxProcessConstraintsCsolnp(&solIneqLB, &solIneqUB, &solEqB);
337
338         if (verbose == 2) {
339             printf("solIneqLB is: ");
340             print(solIneqLB); putchar('\n');
341             printf("solIneqUB is: ");
342             print(solIneqUB); putchar('\n');
343             printf("solEqB is: ");
344             print(solEqB); putchar('\n');
345         }
346     }
347     omxSetupBoundsAndConstraints(freeVarGroup, bl, bu);
348     
349     Matrix blvar = fillMatrix(n, 1, bl);
350     Matrix buvar = fillMatrix(n, 1, bu);
351         
352     /* Initialize Starting Values */
353     if(OMX_VERBOSE) {
354         printf("--------------------------");
355         printf("Starting Values (%d) are:", n);
356     }
357     for(k = 0; k < n; k++) {
358         if((M(myPars, k, 0) == 0.0)) {
359             M(myPars, k, 0) += 0.1;
360         }
361         if(OMX_VERBOSE) { printf("%d: %f", k, M(myPars, k, 0)); }
362     }
363     if(OMX_DEBUG) {
364         printf("--------------------------");
365         printf("Setting up optimizer...");
366     }
367         
368     
369     p_obj = solnp(myPars, solFun, solEqB, solEqBFun, solIneqFun, blvar, buvar, solIneqUB, solIneqLB, myControl, myDEBUG, verbose);
370     
371     
372     fc->fit = p_obj.objValue;
373     if (verbose >= 1) {
374             printf("final objective value is: \n");
375             printf("%2f", fc->fit); putchar('\n');
376     }
377     param_hess = p_obj.parameter;
378     
379     int i;
380     myPars = subset(param_hess, 0, 0, n-1);
381     
382     if (verbose>= 1){
383             printf("final myPars value is: \n");
384             print(myPars); putchar('\n');
385     }
386     myhess = subset(param_hess, 0, n, param_hess.cols - myPars.cols - 2);
387     
388     
389     Matrix inform_m = subset(param_hess, 0, param_hess.cols-1, param_hess.cols-1);
390     
391     inform = M(inform_m, 0, 0);
392     
393     if (verbose >= 2){
394             printf("myhess is: \n");
395             print(myhess); putchar('\n');
396     }
397     
398     double Hess[myhess.cols];
399     
400     for (i = 0; i < myhess.cols; i++)
401     {
402         Hess[i] = myhess.t[i];
403       
404     }
405     
406     mygrad = subset(param_hess, 0, myPars.cols + (myPars.cols*myPars.cols), param_hess.cols-2);
407     
408     
409     memcpy(fc->hess, Hess, sizeof(double) * n * n);
410
411     for (i = 0; i < myPars.cols; i++){
412         x[i] = myPars.t[i];
413     }
414     
415     omxSaveCheckpoint(x, fc->fit, TRUE);
416     
417     GLOB_fc->copyParamToModel(globalState);
418     
419     *inform_out = inform;
420     *iter_out   = iter;
421     
422     GLOB_fitMatrix = NULL;
423     GLOB_fc = NULL;
424     freeMatrices();
425 }
426
427
428
429 void omxCSOLNPConfidenceIntervals(omxMatrix *fitMatrix, FitContext *fc, int verbose)
430 {
431         int ciMaxIterations = Global->ciMaxIterations;
432         // Will fail if we re-enter after an exception
433         //if (NPSOL_fitMatrix) error("NPSOL is not reentrant");
434     
435     GLOB_fitMatrix = fitMatrix;
436         GLOB_fc = fc;
437     
438         FreeVarGroup *freeVarGroup = fitMatrix->fitFunction->freeVarGroup;
439     
440     double inform;
441     double *bl=NULL, *bu=NULL;
442     
443     int n = int(freeVarGroup->vars.size());
444     int ncnln = globalState->ncnln;
445     
446     double optimum = fc->fit;
447     
448     double *optimalValues = fc->est;
449     
450     double f = optimum;
451     std::vector< double > x(n, *optimalValues);
452     std::vector< double > gradient(n);
453     std::vector< double > hessian(n * n);
454     
455     /* CSOLNP Arguments */
456     double EMPTY = -999999.0;
457     
458     Param_Obj p_obj_conf;
459     Matrix param_hess;
460     Matrix myhess = fill(n*n, 1, (double)0.0);
461     Matrix mygrad;
462     Matrix solIneqLB;
463     Matrix solIneqUB;
464     Matrix solEqB;
465     
466     Matrix myPars = fillMatrix(n, 1, fc->est);
467     double (*solFun)(struct Matrix myPars, int verbose);
468     solFun = &csolnpLimitObjectiveFunction;
469     Matrix (*solEqBFun)(int verbose);
470     solEqBFun = &csolnpEqualityFunction;
471     Matrix (*solIneqFun)(int verbose);
472     solIneqFun = &csolnpIneqFun;
473     
474     
475     /* Set boundaries and widths. */
476     /* Allocate arrays */
477     bl      = Calloc ( n + ncnln, double );
478     bu      = Calloc (n + ncnln, double );
479     
480     
481     struct Matrix myControl = fill(6,1,(double)0.0);
482     M(myControl,0,0) = 1.0;
483     M(myControl,1,0) = 400.0;
484     M(myControl,2,0) = 800.0;
485     M(myControl,3,0) = 1.0e-7;
486     M(myControl,4,0) = 1.0e-16;
487     M(myControl,5,0) = 0.0;
488     
489     bool myDEBUG = false;
490     /* Set up actual run */
491     
492     /* needs treatment*/
493     if (ncnln == 0)
494     {
495         solIneqLB = fill(1, 1, EMPTY);
496         solIneqUB = fill(1, 1, EMPTY);
497         solEqB = fill(1, 1, EMPTY);
498     }
499     else{
500         int j;
501         int nineqn;
502         int eqn = 0;
503         for(j = 0; j < globalState->numConstraints; j++) {
504             if (globalState->conList[j].opCode == 1)
505             {
506                 eqn += globalState->conList[j].size;
507             }
508         }
509         if (eqn == ncnln) nineqn = 1;
510         else nineqn = ncnln - eqn;
511         
512         solIneqLB = fill(nineqn, 1, EMPTY);
513         solIneqUB = fill(nineqn, 1, EMPTY);
514             if (eqn == 0) {
515                     solEqB = fill(1, 1, EMPTY);
516             } else {
517                     solEqB = fill(eqn, 1, EMPTY);
518             }
519         
520         omxProcessConstraintsCsolnp(&solIneqLB, &solIneqUB, &solEqB);
521         if (verbose == 2) {
522             printf("solIneqLB is: ");
523             print(solIneqLB); putchar('\n');
524             printf("solIneqUB is: ");
525             print(solIneqUB); putchar('\n');
526             printf("solEqB is: ");
527             print(solEqB); putchar('\n');
528         }
529     }
530     
531     omxSetupBoundsAndConstraints(freeVarGroup, bl, bu);
532     Matrix blvar = fillMatrix(n, 1, bl);
533     Matrix buvar = fillMatrix(n, 1, bu);
534     
535     if(OMX_DEBUG) { mxLog("Calculating likelihood-based confidence intervals."); }
536     
537     
538     if(OMX_DEBUG) { mxLog("Calculating likelihood-based confidence intervals.");
539         mxLog("numIntervals is: %d", Global->numIntervals);
540     }
541     
542     int count = 0;
543     for(int i = 0; i < Global->numIntervals; i++) {
544         if (i == 0)
545         {count++;}
546         omxConfidenceInterval *currentCI = &(Global->intervalList[i]);
547         
548         int msgLength = 45;
549         
550         if (currentCI->matrix->name == NULL) {
551             msgLength += strlen(anonMatrix);
552         } else {
553             msgLength += strlen(currentCI->matrix->name);
554         }
555         
556         char *message = Calloc(msgLength, char);
557         
558         if (currentCI->matrix->name == NULL) {
559             snprintf(message, msgLength, "%s[%d, %d] begin lower interval",
560                      anonMatrix, currentCI->row + 1, currentCI->col + 1);
561         } else {
562             snprintf(message, msgLength, "%s[%d, %d] begin lower interval",
563                      currentCI->matrix->name, currentCI->row + 1, currentCI->col + 1);
564         }
565         
566         omxWriteCheckpointMessage(message);
567         
568         memcpy(x.data(), optimalValues, n * sizeof(double)); // Reset to previous optimum
569         myPars = fillMatrix(n, 1, x.data());
570         CSOLNP_currentInterval = i;
571         
572         
573         currentCI->lbound += optimum;          // Convert from offsets to targets
574         currentCI->ubound += optimum;          // Convert from offsets to targets
575         
576         /* Set up for the lower bound */
577         inform = -1;
578         // Number of times to keep trying.
579         
580         
581         int cycles = ciMaxIterations;
582         
583         double value = INF;
584         double objDiff = 1.e-4;     // TODO : Use function precision to determine CI jitter?
585         
586         while(inform!= 0 && cycles > 0) {
587             /* Find lower limit */
588             currentCI->calcLower = TRUE;
589             
590             p_obj_conf = solnp(myPars, solFun, solEqB, solEqBFun, solIneqFun, blvar, buvar, solIneqUB, solIneqLB, myControl, myDEBUG, verbose);
591             
592             f = p_obj_conf.objValue;
593             
594             myPars = subset(p_obj_conf.parameter, 0, 0, n-1);
595             myhess = subset(p_obj_conf.parameter, 0, n, p_obj_conf.parameter.cols - myPars.cols - 2);
596             
597             mygrad = subset(p_obj_conf.parameter, 0, myPars.cols + (myPars.cols*myPars.cols), p_obj_conf.parameter.cols-2);
598             
599             
600             Matrix inform_m = subset(p_obj_conf.parameter, 0, p_obj_conf.parameter.cols-1, p_obj_conf.parameter.cols-1);
601             
602             inform = M(inform_m, 1, 0);
603             
604             if(verbose>=1) { mxLog("inform_lower is: %f", inform);}
605             currentCI->lCode = inform;
606             
607             
608             if(f < value) {
609                 currentCI->min = omxMatrixElement(currentCI->matrix, currentCI->row, currentCI->col);
610                 
611                 value = f;
612                                 for (int ii = 0; ii < myPars.cols; ii++){
613                                 x.data()[ii] = myPars.t[ii];
614                         }
615                 omxSaveCheckpoint(x.data(), f, TRUE);
616             }
617             
618             if(inform != 0 && OMX_DEBUG) {
619                 mxLog("Calculation of lower interval %d failed: Bad inform value of %f",
620                       i, inform);
621             }
622             cycles--;
623             if(inform != 0) {
624                 unsigned int jitter = TRUE;
625                 for(int j = 0; j < n; j++) {
626                     if(fabs(x[j] - optimalValues[j]) > objDiff) {
627                         jitter = FALSE;
628                         if(OMX_DEBUG) {mxLog("are u here?????");}
629                         break;
630                     }
631                 }
632                 if(jitter) {
633                     for(int j = 0; j < n; j++) {
634                         x[j] = optimalValues[j] + objDiff;
635                     }
636                 }
637             }
638         }
639         
640         if(OMX_DEBUG) { mxLog("Found lower bound %d.  Seeking upper.", i); }
641         // TODO: Repopulate original optimizer state in between CI calculations
642         
643         if (currentCI->matrix->name == NULL) {
644             snprintf(message, msgLength, "%s[%d, %d] begin upper interval",
645                      anonMatrix, currentCI->row + 1, currentCI->col + 1);
646         } else {
647             snprintf(message, msgLength, "%s[%d, %d] begin upper interval",
648                      currentCI->matrix->name, currentCI->row + 1, currentCI->col + 1);
649         }
650         
651         omxWriteCheckpointMessage(message);
652         
653         Free(message);
654         
655         memcpy(x.data(), optimalValues, n * sizeof(double));
656         myPars = fillMatrix(n, 1, x.data());
657         
658         
659         /* Reset for the upper bound */
660         value = INF;
661         inform = -1;
662         cycles = ciMaxIterations;
663         if(verbose >= 1) { mxLog("cycles_upper is: %d", cycles); }
664         while(inform != 0 && cycles >= 0) {
665             /* Find upper limit */
666             currentCI->calcLower = FALSE;
667             p_obj_conf = solnp(myPars, solFun, solEqB, solEqBFun, solIneqFun, blvar, buvar, solIneqUB, solIneqLB, myControl, myDEBUG, verbose);
668             
669             f = p_obj_conf.objValue;
670             
671             myPars = subset(p_obj_conf.parameter, 0, 0, n-1);
672             myhess = subset(p_obj_conf.parameter, 0, n, p_obj_conf.parameter.cols - myPars.cols - 2);
673             
674             mygrad = subset(p_obj_conf.parameter, 0, myPars.cols + (myPars.cols*myPars.cols), p_obj_conf.parameter.cols-2);
675             
676             Matrix inform_m = subset(p_obj_conf.parameter, 0, p_obj_conf.parameter.cols-1, p_obj_conf.parameter.cols-1);
677             
678             
679             inform = M(inform_m, 1, 0);
680             if(verbose >= 1) { mxLog("inform_upper is: %f", inform);}
681             currentCI->uCode = inform;
682             
683             if(f < value) {
684                 currentCI->max = omxMatrixElement(currentCI->matrix, currentCI->row, currentCI->col);
685                 
686                 value = f;
687                                 for (int ii = 0; ii < myPars.cols; ii++){
688                                 x.data()[ii] = myPars.t[ii];
689                 }
690                 
691                 omxSaveCheckpoint(x.data(), f, TRUE);
692             }
693             
694             if(inform != 0 && OMX_DEBUG) {
695                 mxLog("Calculation of upper interval %d failed: Bad inform value of %f",
696                       i, inform);
697             }
698             cycles--;
699             if(inform != 0) {
700                 unsigned int jitter = TRUE;
701                 for(int j = 0; j < n; j++) {
702                     if(fabs(x[j] - optimalValues[j]) > objDiff){
703                         jitter = FALSE;
704                         break;
705                     }
706                 }
707                 if(jitter) {
708                     for(int j = 0; j < n; j++) {
709                         x[j] = optimalValues[j] + objDiff;
710                     }
711                 }
712             }
713         }
714         
715         if(OMX_DEBUG) {mxLog("Found Upper bound %d.", i);}
716         
717     }
718     
719         GLOB_fc = NULL;
720         GLOB_fitMatrix = NULL;
721         CSOLNP_currentInterval = -1;
722 }