Doxygen/OS/_o_s_algorithmic_diff_test_8cpp_source.html

 /* $Id: OSAlgorithmicDiffTest.cpp 2698 2009-06-09 04:14:07Z kmartin $ */

 #include <cstddef>

 #include <cstdlib>

 #include <cctype>

 #include <cassert>

 #include <stack>

 #include <cppad/cppad.hpp>

 #include <iostream>


 // CoinHelperFunctions <has cmath>

 #include "CoinHelperFunctions.hpp"

 #include "OSInstance.h"

 #include "OSiLWriter.h"

 #include "OSParameters.h"

 #include "OSnLNode.h"

 #include "OSErrorClass.h"

 #include "OSFileUtil.h"

 #include "OSiLReader.h"

 #include "OSInstance.h"

 #include "OSExpressionTree.h"

 #include "OSnLNode.h"

 #include "OSDataStructures.h"


 #include <vector>

 #include <map>

 #include <string>


 int  main(){

         WindowsErrorPopupBlocker();

         using std::cout;

         using std::endl;

         using CppAD::AD;

         //using CppAD::NearEqual;

         using CppAD::vector;

         std::cout.precision(12);

         // error checking functions

         bool CheckFunctionValues( double *conVals, double objValue,

                 double x0, double x1, double x2, double x3, double y0, double y1, double z );

         bool CheckHessianUpper( SparseHessianMatrix *sparseHessian,

                 double x0, double x1, double x2, double x3, double y0, double y1, double z );

         bool CheckGradientValues( SparseJacobianMatrix *sparseJac, double *objGrad,

                 double x0, double x1, double x2, double x3, double y0, double y1, double z );

         bool ok = true;

         int k, idx;

         //

         // get the problem data

         //

         FileUtil *fileUtil = NULL;

         std::string osilFileName;

         std::string osil;

         // get the input file

         const char dirsep =  CoinFindDirSeparator();

         // Set directory containing mps data files.

         std::string dataDir;

     dataDir = dirsep == '/' ? "../data/" : "..\\data\\";

         //osilFileName =  dataDir + "osilFiles" + dirsep + "HS071_NLP.osil";

         osilFileName =  dataDir   + "osilFiles" + dirsep + "CppADTestLag.osil";

         std::cout  << "osilFileName  =  " << osilFileName << std::endl;

         fileUtil = new FileUtil();

         osil = fileUtil->getFileAsString( &osilFileName[0]);

         //

         // create OSReader and OSInstance objects

         OSiLReader *osilreader = NULL;

         OSInstance *osinstance = NULL;

         // create reader, parse the OSiL, and generate the OSInstance object

         try{

                 // a counter

                 int kjl;

                 osilreader = new OSiLReader();

                 osinstance = osilreader->readOSiL( osil);

                 std::vector<double> funVals(3);

                 std::vector<double> dfunVals(6);

                 double *conVals = NULL;

                 //conVals = new double[ 2];

                 double *objVals = NULL;

                 //objVals = new double[ 1];


                 //

                 // first initialize the nonlinear structures for call backs

                 std::cout << "Initialize Nonlinear Structures" << std::endl;

                 osinstance->initForAlgDiff( );

                 //osinstance->getJacobianSparsityPattern( );


                 std::map<int, int> varIndexMap;

                 std::map<int, int>::iterator posVarIndexMap;

                 varIndexMap = osinstance->getAllNonlinearVariablesIndexMap( );

                 for(posVarIndexMap = varIndexMap.begin(); posVarIndexMap != varIndexMap.end(); ++posVarIndexMap){

                                 std::cout <<  "Variable Index = "   << posVarIndexMap->first  << std::endl ;

                 }

                 std::cout << "Number of nonlinear variables =  " << varIndexMap.size() << std::endl;


                 //


                 // get the number of nonlinear terms


                 int mm = osinstance->getNumberOfNonlinearExpressionTreeModIndexes();


                 int jj;


                 for(jj = 0; jj < mm; jj++){

                         std::cout << osinstance->getNonlinearExpressionTreeModIndexes()[ jj] << std::endl;

                 }


                 std::cout << "Number of unique nonlinear terms =  " <<  mm << std::endl;

                 //return 0;


                 // domain space vector

                 size_t n  = varIndexMap.size(); // three variables

                 // range space vector

                 size_t m = 3; // Lagrangian has an objective and two constraints


                 std::vector<double> x0( n);

                 x0[0] = 1; // the value for variable x0

                 x0[1] = 5; // the value for variable x1

                 x0[2] = 5; // the value for variable x3

                 std::cout << "CALL forward" << std::endl;

                 funVals = osinstance->forwardAD(0, x0);

                 for( kjl = 0; kjl < 3; kjl++){

                         std::cout << "forward 0 " << funVals[ kjl] << std::endl;

                 }

                 // get the third column of the Jacobian from a forward sweep

                 std::vector<double> x1( n);

                 x1[0] = 0;

                 x1[1] = 0;

                 x1[2] = 1;

                 std::cout << "Now get the third column of the Jacobian forwardAD(1, x1)"  << std::endl;

                 funVals = osinstance->forwardAD(1, x1);

                 for( kjl = 0; kjl < 3; kjl++){

                         std::cout << "forward 1 " << funVals[ kjl] << std::endl;

                 }


                 x1[0] = 1;

                 x1[1] = 0;

                 x1[2] = 1;

                 // recalculate the forward sweep with the new x1 vector

                 funVals = osinstance->forwardAD(1, x1);

                 std::vector<double> x2( n);

                 x2[0] = 0;

                 x2[1] = 0;

                 x2[2] = 0;

                 std::cout << "Now calcuate forwardAD(2, x2)"  << std::endl;

                 funVals = osinstance->forwardAD(2, x2);

                 for( kjl = 0; kjl < 3; kjl++){

                         std::cout << "forward 2 " << funVals[ kjl] << std::endl;

                 }


                 std::vector<double> vlambda(3);

                 vlambda[0] = 0;

                 vlambda[1] = 0;

                 vlambda[2] = 1;

                 // reverse sweep to get third row of Jacobian

                 std::cout << "Now get the third row of the Jacobian reverseAD(1, vlambda)"  << std::endl;

                 osinstance->forwardAD(0, x0);

                 funVals = osinstance->reverseAD(1, vlambda);

                 for( kjl = 0; kjl < 3; kjl++){

                         std::cout << "reverse 1 " << funVals[ kjl] << std::endl;

                 }

                 // now get the Hessian of the Lagrangian of objective and

                 // with the following multipliers

                 vlambda[0] = 1;

                 vlambda[1] = 2;

                 vlambda[2] = 1;

                 x1[0] = 0;

                 x1[1] = 0;

                 x1[2] = 1;

                 funVals = osinstance->forwardAD(1, x1);


                 dfunVals = osinstance->reverseAD(2, vlambda);

                 // get the first partials of the Lagrangian

                 std::cout << "Here are the first partials of the Lagrangain" << std::endl;

                 for(int kjl = 0; kjl <= 4; kjl+=2){

                         std::cout << dfunVals[ kjl] << std::endl;

                 }

                 std::cout << "Here is the third row (column) of Hessian of Lagrangian" << std::endl;

                 for(int kjl = 1; kjl <= 5; kjl+=2){

                         std::cout << dfunVals[ kjl] << std::endl;

                 }

                 double* x = new double[4]; //primal variables

                 double* z = new double[2]; //Lagrange multipliers on constraints

                 double* w = new double[1]; //Lagrange multiplier on objective

                 x[ 0] = 1;    // primal variable 0

                 x[ 1] = 5;    // primal variable 1

                 x[ 2] = 10;   // primal variable 2

                 x[ 3] = 5;    // primal variable 3

                 z[ 0] = 2;    // Lagrange multiplier on constraint 0

                 z[ 1] = 1;    // Lagrange multiplier on constraint 1

                 w[ 0] = 1;    // Lagrange multiplier on the objective function


                 osinstance->bUseExpTreeForFunEval = false;

                 std::cout << "Calculate objective, idx = -1"  << std::endl;

                 std::cout << "obj value = " << osinstance->calculateFunctionValue(-1, x,  true) << std::endl;


                 std::cout << "Calculate  first constraint, idx = 0"  << std::endl;

                 std::cout << "constraint index 0 value = " << osinstance->calculateFunctionValue(0, x,  true) << std::endl;


                 std::cout << "Now use calculateAllConstraintFunctionValues"  << std::endl;

                 conVals = osinstance->calculateAllConstraintFunctionValues(x, true);

                 for( idx = 0; idx < osinstance->getConstraintNumber(); idx++){

                         std::cout << "CONSTRAINT FUNCTION INDEX = " <<  idx << "  CONSTRAINT FUNCTION VALUE =  "  << *(conVals + idx) << std::endl;

                 }

                 //

                 std::cout << "Now use calculateAllObjectiveFunctionValues"  << std::endl;

                 objVals = osinstance->calculateAllObjectiveFunctionValues( x, NULL, NULL, true, 0);

                 for( idx = 0; idx < osinstance->getObjectiveNumber(); idx++){

                         std::cout << "OBJECTIVE FUNCTION  INDEX = " << idx <<  "  OBJECTIVE FUNCTION VALUE = "  << *(objVals + idx) << std::endl;

                 }

                 ok = CheckFunctionValues( conVals, *objVals, x[ 0], x[1], x[2], x[3],  z[0], z[1], w[0] );

                 if( ok == 0){

                         std::cout << "FAILED CHECKING FUNCTION VALUES TEST" << std::endl;

                         return 0;

                 }

                 else{

                         std::cout << "PASSED CHECKING FUNCTION VALUES TEST" << std::endl;

                 }


                 std::cout << "PLACE CALL TO JACOBIAN SPARSITY PATTERN"   << std::endl;

                 SparseJacobianMatrix *sparseJac;

                 sparseJac = osinstance->getJacobianSparsityPattern();

                 // print out just the sparsity pattern

                 std::cout << "JACOBIAN SPARSITY PATTERN"   << std::endl;

                 std::cout << "JACOBIAN START SIZE "   <<  sparseJac->startSize << std::endl;

                 for(idx = 0; idx < osinstance->getConstraintNumber(); idx++){

                         // some solvers (e.g. reduced gradient solvers) may want to know which values of the

                         // Jacobian matrix are constant, i.e. linear, sparseJac->conVals is the number of constant

                         // terms in the gradient for each rowt, the first conVals terms are constant, when getting

                         std::cout << "number constant terms in constraint "   <<  idx << " is "

                         << *(sparseJac->conVals + idx)  << std::endl;

                         for(k = *(sparseJac->starts + idx); k < *(sparseJac->starts + idx + 1); k++){

                                 std::cout << "row idx = " << idx <<  "  col idx = "<< *(sparseJac->indexes + k) << std::endl;

                         }

                 }


                 SparseHessianMatrix *sparseHessian;

                 // the Hessian test

                 // get the sparsity pattern -- many solvers want to initialize with just the sparsity

                 std::cout << "GET LAGRANGIAN HESSIAN SPARSITY PATTERN"   << std::endl;

                 sparseHessian = osinstance->getLagrangianHessianSparsityPattern( );

                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){

                         std::cout <<  "Row Index = " << *(sparseHessian->hessRowIdx + idx) ;

                         std::cout <<  "  Column Index = " << *(sparseHessian->hessColIdx + idx) << std::endl;

                 }


                 double *objGrad;

                 std::cout << "OBJECTIVE FUNCTION GRADIENT"   << std::endl;

                 // in our implementation the objective function is a dense gradient

                 objGrad = osinstance->calculateObjectiveFunctionGradient( x, NULL, NULL,  -1, false, 1);

                 for(idx = 0; idx < osinstance->getVariableNumber(); idx++){

                         std::cout << "col idxx = " << idx << "  value =  " << *(objGrad + idx)  << std::endl;

                 }

                 std::cout << "CONSTRAINT JACOBIAN MATRIX"   << std::endl;

                 // now make the gradient calculations and fill in the sparse Jacobian matrix

                 sparseJac = osinstance->calculateAllConstraintFunctionGradients( x, NULL, NULL,  false, 1);

                 for(idx = 0; idx < osinstance->getConstraintNumber(); idx++){

                         for(k = *(sparseJac->starts + idx); k < *(sparseJac->starts + idx + 1); k++){

                                 std::cout << "row idx = " << idx <<  "  col idx = "<< *(sparseJac->indexes + k)

                                 << " value = " << *(sparseJac->values + k) << std::endl;

                         }

                 }

                 ok = CheckGradientValues( sparseJac, objGrad, x[ 0], x[1], x[2], x[3], z[0], z[1], w[0] );

                 if( ok == 0){

                         std::cout << "FAILED THE GRADIENT TEST" << std::endl;

                         return 0;

                 }

                 else{

                         std::cout << "PASSED THE GRADIENT TEST" << std::endl;

                 }

                 //first iteration

                 std::cout << "GET LAGRANGIAN HESSIAN FIRST TIME"   << std::endl;

                 sparseHessian = osinstance->calculateLagrangianHessian( x, w,  z,  false, 2);

                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){

                         std::cout << "row idx = " << *(sparseHessian->hessRowIdx + idx) <<

                         "  col idx = "<< *(sparseHessian->hessColIdx + idx)

                         << " value = " << *(sparseHessian->hessValues + idx) << std::endl;

                 }

                 ok = CheckHessianUpper( sparseHessian, x[0],  x[1], x[2],  x[3],  z[0], z[1], w[0]);

                 if( ok == 0){

                         std::cout << "FAILED THE FIRST HESSIAN TEST" << std::endl;

                         return 0;

                 }

                 else{

                         std::cout << "PASSED THE FIRST HESSIAN TEST" << std::endl;

                 }

                 // now change an x value, we don't rebuild the tree, however new_x

                 // must be set to true

                 x[0] = 5;


                 std::cout << "NOW GET LAGRANGIAN HESSIAN SECOND TIME FOR x[0] = 5"   << std::endl;


                 osinstance->calculateAllObjectiveFunctionValues( x, w, z, true, 2);

                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){

                         std::cout << "row idx = " << *(sparseHessian->hessRowIdx + idx) <<

                         "  col idx = "<< *(sparseHessian->hessColIdx + idx)

                         << " value = " << *(sparseHessian->hessValues + idx) << std::endl;

                 }

                 ok = CheckHessianUpper( sparseHessian , x[0],  x[1], x[2], x[3],  z[0], z[1], w[0] );

                 if( ok == 0){

                         std::cout << "FAILED THE SECOND HESSIAN TEST" << std::endl;

                         return 0;

                 }

                 else{

                         std::cout << "PASSED THE SECOND HESSIAN TEST" << std::endl  << std::endl ;

                 }

                 std::cout << "HERE IS ROW 1 OF JACOBIAN MATRIX" << std::endl;

                 idx = 1;

                 for(k = *(sparseJac->starts + idx); k < *(sparseJac->starts + idx + 1); k++){

                         std::cout << "row idx = " << idx <<  "  col idx = "<< *(sparseJac->indexes + k)

                                 << " value = " << *(sparseJac->values + k) << std::endl;

                 }

                 std::cout << std::endl;

                 sparseHessian = osinstance->calculateHessian(x, 1, true);

                 std::cout << "HERE IS ROW 1 HESSIAN MATRIX" << std::endl;

                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){

                         std::cout << "row idx = " << *(sparseHessian->hessRowIdx + idx) <<

                         "  col idx = "<< *(sparseHessian->hessColIdx + idx)

                         << " value = " << *(sparseHessian->hessValues + idx) << std::endl;

                 }

                 //

                 // adjust the Lagrange multipliers to correspond to finding Hessian of constraint 1

                 z[ 0] = 0;  // Lagrange multiplier on constraint 0

                 z[ 1] = 1;  // Lagrange multiplier on constraint 1

                 w[ 0] = 0;  // Lagrange multiplier on the objective function

                 ok = CheckHessianUpper( sparseHessian , x[0],  x[1], x[2], x[3], z[0], z[1], w[0] );

                 if( ok == 0){

                         std::cout << "FAILED THE THIRD HESSIAN TEST" << std::endl;

                         return 0;

                 }

                 else{

                         std::cout << "PASSED THE THIRD HESSIAN TEST" << std::endl  << std::endl ;

                 }

                 //set x[0] back to its original value of 1

                 x[ 0] = 1;

                 //return 0;

                 //

                 //

                 // now work directly with the CppAD package instead of OSInstance API

                 //

                 n = 4;

                 m = 3;

                 CppADvector< AD<double> > X(n);

                 CppADvector< AD<double> > Y(m);

                 X[0] = 5;

                 X[1] = 5;

                 X[2] = 0;

                 X[3] = 1;

                 // declare independent variables and start tape recording

                 std::cout << "Start Taping" << std::endl;

                 CppAD::Independent( X);

                 // range space vector

                 // we include the constant terms in the CppAD functions

                 Y[ 0] =  CppAD::pow(X[0], 2) + 9*X[1];

                 Y[ 1] =  33 - 105 + 1.37*X[1] + 2*X[3] + 5*X[1] ;

                 Y[ 2] =  log(X[0]*X[3]) + 7*X[2] ;

                 // create f: x -> y and stop tape recording

                 CppAD::ADFun<double> f(X, Y);

                 std::cout << "Stop Taping" << std::endl;

                 // get function values

                 std::vector<double> x_vec( n);

                 x_vec[ 0] = x[ 0];

                 x_vec[ 1] = x[ 1];

                 x_vec[ 2] = x[ 2];

                 x_vec[ 3] = x[ 3];

                 funVals = f.Forward(0, x_vec);

                 conVals[ 0] = funVals[ 1];

                 std::cout << "conVals[ 0] = " << conVals[ 0] << std::endl;

                 conVals[ 1] = funVals[ 2];

                 std::cout << "conVals[ 1] = " << conVals[ 1] << std::endl;

                 objVals[ 0] = funVals[ 0];

                 std::cout << "objVals[ 0] = " << objVals[ 0] << std::endl;

                 ok = CheckFunctionValues( conVals, funVals[ 0], x[ 0], x[1], x[2], x[3], z[0], z[1], w[0] );

                 if( ok == 0){

                         std::cout << "FAILED CHECKING FUNCTION VALUES TEST" << std::endl;

                         return 0;

                 }

                 else{

                         std::cout << "PASSED CHECKING FUNCTION VALUES TEST" << std::endl;

                 }

                 // now get gradient and Hessian

                 // first define and initialze unit vector vector

                 sparseJac = osinstance->getJacobianSparsityPattern();

                 std::vector<double> unit_col_vec( n);

                 std::vector<double> lagMultipliers( m);

                 std::vector<double> gradVals( m);

                 lagMultipliers[ 0] = w[ 0];

                 lagMultipliers[ 1] = z[ 0];

                 lagMultipliers[ 2] = z[ 1];

                 unsigned int index, kj;

                 //return 0;

                 for(index = 0; index < n; index++){

                         unit_col_vec[ index] = 0;

                 }

                 for(index = 0; index < n; index++){

                         unit_col_vec[ index] = 1;

                         // calculate column i of the Jacobian matrix

                         gradVals = f.Forward(1, unit_col_vec);

                         unit_col_vec[ index] = 0;

                         // get the nonzero gradient values in constraint k

                         for(kj = 0; kj < m; kj++){

                                 std::cout << "variable " << index << "  row " << kj << "  gradient value" << std::endl;

                                 std::cout << "gradient value = " << gradVals[ kj] << std::endl;

                         }

                         // get row i of the Lagrangian function!!!

                         std::cout << "CALL f.Reverse -------" << std::endl;

                         f.Reverse(2, lagMultipliers);

                         std::cout << "FINISH CALL f.Reverse -------" << std::endl;

                 }

                 // done with CppAD test

                 // do garbage collection

                 delete osilreader;

                 osilreader = NULL;

                 std::cout << "OSILREADER DELETED" << std::endl;

                 //delete[] conVals;

                 //delete[] objVals;

                 delete[] x;

                 delete[] z;

                 delete[] w;

         }

         catch(const ErrorClass& eclass){

                 std::cout << eclass.errormsg << std::endl;

         }


         //

         {

                 //checking CppAD power, another illustration of CppAD

                 size_t n  = 2;

                 double x0 = 4;

                 double x1 = .5;

                 CppADvector< AD<double> > x(n);

             x[0]      = x0;

             x[1]      = x1;

              // declare independent variables and start tape recording

              CppAD::Independent(x);

              // range space vector

              size_t m = 1;

              CppADvector< AD<double> > y(m);

              y[0] = std::pow(x0, x1);

              // create f: x -> y and stop tape recording

              CppAD::ADFun<double> f(x, y);

              // check value

              double check = std::pow(x0, x1);

              // forward computation of first partial w.r.t. x[0]

              std::vector<double> dx(n);

              std::vector<double> dy(m);

              dx[0] = 4.;

              dx[1] = 1/2.;

              dy    = f.Forward(1, dx);

              std::cout << "dy =  " <<  dy[ 0] << std::endl;

              check = x1 * std::pow(x0, x1-1.);

              //ok   &= NearEqual(dy[0], check, 1e-10, 1e-10);

              ok = ( fabs(check - dy[0])/(fabs( check) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         }


         {


                 //checking CppAD sparsity features

                 // domain space vector

                 size_t n = 3;

                 CPPAD_TEST_VECTOR< AD<double> > X(n);

                 X[0] = 0.;

                 X[1] = 1.;

                 X[2] = 2.;

                 // declare independent variables and start recording

                 CppAD::Independent(X);

                 // range space vector

                 size_t m = 2;

                 CPPAD_TEST_VECTOR< AD<double> > Y(m);

                 Y[0] = CppAD::pow(X[0], 2)  + CppAD::pow(X[2], 2);

                 Y[1] = -CppAD::pow(X[0], 2) + CppAD::pow(X[1], 2);

                 // create f: X -> Y and stop tape recording

                 CppAD::ADFun<double> f(X, Y);


                 // sparsity pattern for the identity matrix

                 std::vector<bool> r(n * n);

                 size_t i, j;

                 for(i = 0; i < n; i++) {

                         for(j = 0; j < n; j++)

                                 r[ i * n + j ] = false;

                                 r[ i * n + i ] = true;

                 }

                 // compute sparsity pattern for J(x) = F^{(1)} (x)

                 f.ForSparseJac(n, r);

                 //std::vector<bool> s(m * m);

                 //for(i = 0; i < m; i++){

                 //      for(j = 0; j < m; j++)

                 //              s[ i * m + j ] = false;

                 //      s[ i * m + i ] = true;

             // }

             // sparsity pattern for F'(x)

             // f.RevSparseJac(m, s);

                 // compute sparsity pattern for H(x) = F_0^{(2)} (x)

                 std::vector<bool> e( m);

                 //Vector s(m);

                 for(i = 0; i < m; i++)

                 e[i] = false;

                 e[ 0] = true;

                 e[ 1] = false;

                 std::vector<bool> h( n*n);

                 //Vector h(n * n);

                 std::cout << "Computing Sparse Hessian" << std::endl;

                 h = f.RevSparseHes(n, e);

                 for(i = 0; i < n; i++){

                         std::cout << "Row " << i << "  of Hessian " << std::endl;

                         for(j = 0; j < n; j++){

                                 std::cout << h[ i*n + j] <<  "  ";

                         }

                         std::cout << std::endl;

                 }

         }

         delete fileUtil;

         std::cout << "\nTEST OF ALGORITHMIC DIFFERENTIATION CONCLUDED SUCCESSFULLY\n";

         return 0;

 }// end main program


 bool CheckFunctionValues( double *conVals, double objValue,

         double x0, double x1, double x2, double x3, double z0, double z1, double w ){

         using CppAD::NearEqual;

         bool ok  = true;

         double checkObj = x0*x0 + 9*x1;

         std::cout  << "checkObj = " << checkObj << std::endl;

         std::cout  << "objValue = " << objValue << std::endl;

         //ok &= NearEqual(objValue, checkObj, 1e-10, 1e-10);

         ok = ( fabs(checkObj - objValue )/(fabs( checkObj) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         double checkCon0 = 33. - 105. + 1.37*x1 + 2*x3 + 5*x1;

         std::cout  << "checkCon0 = " << checkCon0 << std::endl;

         std::cout  << "conVals = " << *(conVals + 0) << std::endl;

         //ok &= NearEqual(*(conVals + 0), checkCon0, 1e-10, 1e-10);

         ok = ( fabs(checkCon0 - *(conVals + 0) )/(fabs( checkCon0) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         double checkCon1 = log(x0*x3) + 7*x2;

         std::cout  << "checkCon1 = " << checkCon1 << std::endl;

         std::cout  << "conVals = " << *(conVals + 1) << std::endl;

         //ok &= NearEqual( *(conVals + 1), checkCon1, 1e-10, 1e-10);

         ok = ( fabs(checkCon1 - *(conVals + 1) )/(fabs( checkCon1) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         return ok;

 }//CheckFunctionValues

 //

 //

 bool CheckGradientValues( SparseJacobianMatrix *sparseJac, double *objGrad,

         double x0, double x1, double x2, double x3, double y0, double y1, double w ){

         using CppAD::NearEqual;

         bool ok  = true;

         // first the objective function gradient

         double checkObjPartial0 = 2*x0;

         //ok &= NearEqual( *(objGrad + 0), checkObjPartial0, 1e-10, 1e-10);

         ok = ( fabs(checkObjPartial0 - *(objGrad + 0) )/(fabs( checkObjPartial0) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         double checkObjPartial1 = 9;

         //ok &= NearEqual( *(objGrad + 1), checkObjPartial1, 1e-10, 1e-10);

         ok = ( fabs(checkObjPartial1 - *(objGrad + 1) )/(fabs( checkObjPartial1) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         double checkObjPartial2 = 0;

         //ok &= NearEqual( *(objGrad + 2), checkObjPartial2, 1e-10, 1e-10);

         ok = ( fabs(checkObjPartial2 - *(objGrad + 2) )/(fabs( checkObjPartial2) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         // get the constrating gradient

         // row 0 gradient -- there are nonzero partials for variables 1 and 2

         double checkCon0Partial1 = 1.37 + 5.0;

         //ok &= NearEqual( *(sparseJac->values + 0), checkCon0Partial1, 1e-10, 1e-10);

         ok = ( fabs(checkCon0Partial1 - *(sparseJac->values + 0) )/(fabs( checkCon0Partial1) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         double checkCon0Partial3 = 2.;

         //ok &= NearEqual( *(sparseJac->values + 1), checkCon0Partial3, 1e-10, 1e-10);

         ok = ( fabs(checkCon0Partial3 - *(sparseJac->values + 1) )/(fabs( checkCon0Partial3) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         // row 1 gradient -- there are nonzero partials for variables 0 and 2

         double checkCon1Partial2 = 7;

         //ok &= NearEqual( *(sparseJac->values + 2), checkCon1Partial2, 1e-10, 1e-10);

         ok = ( fabs(checkCon1Partial2 - *(sparseJac->values + 2) )/(fabs( checkCon1Partial2) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         double checkCon1Partial0 = 1./x0;

         //ok &= NearEqual( *(sparseJac->values + 3), checkCon1Partial0, 1e-10, 1e-10);

         ok = ( fabs(checkCon1Partial0 - *(sparseJac->values + 3) )/(fabs( checkCon1Partial0) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         double checkCon1Partial3 = 1./x3;

         //ok &= NearEqual( *(sparseJac->values + 4), checkCon1Partial3, 1e-10, 1e-10);

         ok = ( fabs(checkCon1Partial3 - *(sparseJac->values + 4) )/(fabs( checkCon1Partial3) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;

         return ok;

 }//CheckGradientValues

 //

 bool CheckHessianUpper( SparseHessianMatrix *sparseHessian ,

         double x0, double x1, double x2, double x3, double z0, double z1, double w ){

         using CppAD::NearEqual;

         bool ok  = true;

         int hessValuesIdx = 0;

         //assert( sparseHessian->hessDimension = n * (n + 1) /2)

         // L_00 = 2 * w - z1 / ( x0 * x0 )

         double check = 2. * w - z1 / (x0 * x0);

         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), check, 1e-10, 1e-10);

         if(ok == false) std::cout << "FAILED ONE" << std::endl;

         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), 0., 1e-10, 1e-10);

         if(ok == false) std::cout << "FAILED TWO" << std::endl;

         // L_22 = - z1 / (x3 * x3)

         check = - z1 / (x3 * x3);

         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), check, 1e-10, 1e-1);

         if(ok == false) std::cout << "FAILED THREE" << std::endl;

         return ok;

 }//CheckHessianUpper

OSErrorClass.h

CheckFunctionValues
bool CheckFunctionValues(double *conVals, double objValue, double x0, double x1, double x2, double x3, double z0, double z1, double w)
Definition: OSAlgorithmicDiffTest.cpp:753

OSnLNode.h
This file defines the OSnLNode class along with its derived classes.

OSInstance::calculateLagrangianHessian
SparseHessianMatrix * calculateLagrangianHessian(double *x, double *objLambda, double *conLambda, bool new_x, int highestOrder)
Calculate the Hessian of the Lagrangian Expression Tree This method will build the CppAD expression t...
Definition: OSInstance.cpp:6135

SparseJacobianMatrix::values
double * values
values holds a double array of nonzero partial derivatives
Definition: OSGeneral.h:340

OSInstance::calculateFunctionValue
double calculateFunctionValue(int idx, double *x, bool new_x)
Calculate the function value for function (constraint or objective) indexed by idx.
Definition: OSInstance.cpp:5769

OSInstance::bUseExpTreeForFunEval
bool bUseExpTreeForFunEval
bUseExpTreeForFunEval is set to true if you wish to use the OS Expression Tree for function evaluatio...
Definition: OSInstance.h:4899

OSInstance::getVariableNumber
int getVariableNumber()
Get number of variables.
Definition: OSInstance.cpp:2018

x2
ULong x2
Definition: OSdtoa.cpp:1776

OSInstance::calculateHessian
SparseHessianMatrix * calculateHessian(double *x, int idx, bool new_x)
Calculate the Hessian of a constraint or objective function.
Definition: OSInstance.cpp:6153

ErrorClass::errormsg
std::string errormsg
errormsg is the error that is causing the exception to be thrown
Definition: OSErrorClass.h:42

OSExpressionTree.h

main
int main(int argc, char *argv[])
Definition: BB_tm.cpp:32

OSInstance::getNonlinearExpressionTreeModIndexes
int * getNonlinearExpressionTreeModIndexes()
Get all the nonlinear expression tree indexes, i.e., indexes of rows (objectives or constraints) that...
Definition: OSInstance.cpp:2720

SparseJacobianMatrix::startSize
int startSize
startSize is the dimension of the starts array – should equal number of rows + 1
Definition: OSGeneral.h:313

OSInstance::calculateAllConstraintFunctionGradients
SparseJacobianMatrix * calculateAllConstraintFunctionGradients(double *x, double *objLambda, double *conLambda, bool new_x, int highestOrder)
Calculate the gradient of all constraint functions.
Definition: OSInstance.cpp:5912

OSiLReader.h

OSFileUtil.h

OSInstance::getNumberOfNonlinearExpressionTreeModIndexes
int getNumberOfNonlinearExpressionTreeModIndexes()
Get the number of unique nonlinear expression tree indexes after modifying the expression tree to con...
Definition: OSInstance.cpp:2714

SparseHessianMatrix::hessColIdx
int * hessColIdx
hessColIdx is an integer array of column indices in the range 0, ..., n - 1.
Definition: OSGeneral.h:399

OSInstance::getObjectiveNumber
int getObjectiveNumber()
Get number of objectives.
Definition: OSInstance.cpp:2174

CheckHessianUpper
bool CheckHessianUpper(SparseHessianMatrix *sparseHessian, double x0, double x1, double x2, double x3, double z0, double z1, double w)
Definition: OSAlgorithmicDiffTest.cpp:811

OSiLReader::readOSiL
OSInstance * readOSiL(const std::string &osil)
parse the OSiL model instance.
Definition: OSiLReader.cpp:53

OSInstance::getAllNonlinearVariablesIndexMap
std::map< int, int > getAllNonlinearVariablesIndexMap()
Definition: OSInstance.cpp:6595

SparseHessianMatrix::hessRowIdx
int * hessRowIdx
hessRowIdx is an integer array of row indices in the range 0, ..., n - 1.
Definition: OSGeneral.h:394

j
static char * j
Definition: OSdtoa.cpp:3622

f
void fint fint fint real fint real real real real * f
Definition: BonBqpdSolver.cpp:27

OSInstance::getJacobianSparsityPattern
SparseJacobianMatrix * getJacobianSparsityPattern()
Definition: OSInstance.cpp:5058

e
void fint fint fint real fint real real real real real real real real real * e
Definition: BonBqpdSolver.cpp:27

x0
ULong * x0
Definition: OSdtoa.cpp:1776

x1
ULong x1
Definition: OSdtoa.cpp:1776

OSiLReader
Used to read an OSiL string.
Definition: OSiLReader.h:37

SparseJacobianMatrix::indexes
int * indexes
indexes holds an integer array of variable indices.
Definition: OSGeneral.h:335

OSInstance::calculateAllConstraintFunctionValues
double * calculateAllConstraintFunctionValues(double *x, double *objLambda, double *conLambda, bool new_x, int highestOrder)
Calculate all of the constraint function values.
Definition: OSInstance.cpp:5834

SparseHessianMatrix::hessDimension
int hessDimension
hessDimension is the number of nonzeros in each array.
Definition: OSGeneral.h:389

OSInstance::getLagrangianHessianSparsityPattern
SparseHessianMatrix * getLagrangianHessianSparsityPattern()
Definition: OSInstance.cpp:6641

k
void fint fint * k
Definition: BonBqpdSolver.cpp:27

SparseHessianMatrix::hessValues
double * hessValues
hessValues is a double array of the Hessian values.
Definition: OSGeneral.h:404

OSDataStructures.h

r
void fint fint fint real fint real real real real real real real * r
Definition: BonBqpdSolver.cpp:27

SparseJacobianMatrix
a sparse Jacobian matrix data structure
Definition: OSGeneral.h:300

OS_NEAR_EQUAL
#define OS_NEAR_EQUAL
we use OS_NEAR_EQUAL in unitTest to see if we are close to the optimal obj value
Definition: OSParameters.h:89

SparseJacobianMatrix::starts
int * starts
starts holds an integer array of start elements, each start element points to the start of partials f...
Definition: OSGeneral.h:324

OSInstance::getConstraintNumber
int getConstraintNumber()
Get number of constraints.
Definition: OSInstance.cpp:2327

FileUtil::getFileAsString
std::string getFileAsString(const char *fname)
read a file and return contents as a string.
Definition: OSFileUtil.cpp:35

OSInstance::initForAlgDiff
bool initForAlgDiff()
This should be called by nonlinear solvers using callback functions.
Definition: OSInstance.cpp:7107

OSInstance::calculateAllObjectiveFunctionValues
double * calculateAllObjectiveFunctionValues(double *x, double *objLambda, double *conLambda, bool new_x, int highestOrder)
Calculate all of the objective function values.
Definition: OSInstance.cpp:5873

OSInstance::reverseAD
std::vector< double > reverseAD(int p, std::vector< double > vdlambda)
Perform an AD reverse sweep.
Definition: OSInstance.cpp:7742

m
void fint * m
Definition: BonBqpdSolver.cpp:27

OSInstance
The in-memory representation of an OSiL instance..
Definition: OSInstance.h:2262

osinstance
OSInstance * osinstance
Definition: OSParseosil.tab.cpp:3860

SparseJacobianMatrix::conVals
int * conVals
conVals holds an integer array of integers, conVals[i] is the number of constant terms in the gradien...
Definition: OSGeneral.h:330

FileUtil
class used to make it easy to read and write files.
Definition: OSFileUtil.h:37

w
void fint fint fint real fint real real real real real real real real * w
Definition: BonBqpdSolver.cpp:27

OSiLWriter.h

SparseHessianMatrix
The in-memory representation of a SparseHessianMatrix..
Definition: OSGeneral.h:376

CheckGradientValues
bool CheckGradientValues(SparseJacobianMatrix *sparseJac, double *objGrad, double x0, double x1, double x2, double x3, double y0, double y1, double w)
Definition: OSAlgorithmicDiffTest.cpp:776

ErrorClass
used for throwing exceptions.
Definition: OSErrorClass.h:31

n
void fint * n
Definition: BonFilterSolver.cpp:96

OSInstance::forwardAD
std::vector< double > forwardAD(int p, std::vector< double > vdX)
Perform an AD forward sweep.
Definition: OSInstance.cpp:7717

OSInstance::calculateObjectiveFunctionGradient
double * calculateObjectiveFunctionGradient(double *x, double *objLambda, double *conLambda, int objIdx, bool new_x, int highestOrder)
Calculate the gradient of the objective function indexed by objIdx.
Definition: OSInstance.cpp:6017

x
void fint fint fint real fint real * x
Definition: BonBqpdSolver.cpp:27

log
CouExpr & log(CouExpr &e)