/home/coin/SVN-release/OS-1.0.0/OS/examples/algorithmicDiff/algorithmicDiffTest.cpp

Go to the documentation of this file.
00001 
00072 #include <cstddef>
00073 #include <cstdlib>
00074 #include <cctype>
00075 #include <cassert>
00076 #include <stack>
00077 #include <cppad/cppad.hpp>
00078 #include <iostream>
00079 
00080 // CoinHelperFunctions <has cmath>
00081 #include "CoinHelperFunctions.hpp"
00082 #include "OSInstance.h"
00083 #include "OSiLWriter.h"
00084 #include "OSParameters.h"
00085 #include "OSnLNode.h"
00086 #include "ErrorClass.h"
00087 #include "FileUtil.h"  
00088 #include "OSiLReader.h"
00089 #include "OSInstance.h"
00090 #include "OSExpressionTree.h"
00091 #include "OSnLNode.h"
00092 #include "OSDataStructures.h"
00093 
00094 
00095 #include <vector>  
00096 #include <map> 
00097 #include<string>
00098  
00099 
00100 int  main(){    
00101         using std::cout;
00102         using std::endl;
00103         using CppAD::AD;
00104         using CppAD::NearEqual;
00105         using CppAD::vector;
00106         std::cout.precision(12);
00107         // error checking functions
00108         bool CheckFunctionValues( double *conVals, double objValue,
00109                 double x0, double x1, double x2, double x3, double y0, double y1, double z );
00110         bool CheckHessianUpper( SparseHessianMatrix *sparseHessian, 
00111                 double x0, double x1, double x2, double x3, double y0, double y1, double z );
00112         bool CheckGradientValues( SparseJacobianMatrix *sparseJac, double *objGrad,
00113                 double x0, double x1, double x2, double x3, double y0, double y1, double z );
00114         bool ok = true;
00115         int k, idx;
00116         //
00117         // get the problem data
00118         //
00119         FileUtil *fileUtil = NULL; 
00120         std::string osilFileName;
00121         std::string osil;
00122         // get the input file
00123         const char dirsep =  CoinFindDirSeparator();
00124         // Set directory containing mps data files.
00125         std::string dataDir;
00126     dataDir = dirsep == '/' ? "../../data/" : "..\\..\\data\\";
00127         //osilFileName =  dataDir + "HS071_NLP.osil";
00128         osilFileName =  dataDir + "osilFiles" + dirsep  + "CppADTestLag.osil";
00129         std::cout  << "osilFileName  =  " << osilFileName << std::endl;
00130         fileUtil = new FileUtil();
00131         osil = fileUtil->getFileAsString( &osilFileName[0]);    
00132         //
00133         // create OSReader and OSInstance objects
00134         OSiLReader *osilreader = NULL;
00135         OSInstance *osinstance = NULL;
00136         // create reader, parse the OSiL, and generate the OSInstance object
00137         try{    
00138                 // a counter
00139                 int kjl;
00140                 osilreader = new OSiLReader();
00141                 osinstance = osilreader->readOSiL( &osil);
00142                 std::vector<double> funVals(3);
00143                 std::vector<double> dfunVals(6);
00144                 double *conVals = NULL;
00145                 conVals = new double[ 2];
00146                 double *objVals = NULL;
00147                 objVals = new double[ 1];
00148 
00162                 //
00163                 // first initialize the nonlinear structures for call backs
00164                 std::cout << "Initialize Nonlinear Structures" << std::endl;
00165                 osinstance->initForAlgDiff( );
00166                 
00172                 std::map<int, int> varIndexMap;
00173                 std::map<int, int>::iterator posVarIndexMap;
00174                 varIndexMap = osinstance->getAllNonlinearVariablesIndexMap( );
00179                 for(posVarIndexMap = varIndexMap.begin(); posVarIndexMap != varIndexMap.end(); ++posVarIndexMap){
00180                                 std::cout <<  "Variable Index = "   << posVarIndexMap->first  << std::endl ;
00181                 }
00182                 std::cout << "Number of nonlinear variables =  " << varIndexMap.size() << std::endl;
00183                 // domain space vector
00184                 size_t n  = varIndexMap.size(); // three variables
00185                 // range space vector
00186                 size_t m = 3; // Lagrangian has an objective and two constraints
00187 
00188                 std::vector<double> x0( n);
00193                 x0[0] = 1; // the value for variable x0
00194                 x0[1] = 5; // the value for variable x1
00195                 x0[2] = 5; // the value for variable x3                         
00196                 std::cout << "CALL forward" << std::endl;
00197                 funVals = osinstance->forwardAD(0, x0);
00198                 for( kjl = 0; kjl < 3; kjl++){
00199                         std::cout << "forward 0 " << funVals[ kjl] << std::endl;
00200                 }
00201                 // get the third column of the Jacobian from a forward sweep
00202                 std::vector<double> x1( n);
00203                 x1[0] = 0;
00204                 x1[1] = 0;
00205                 x1[2] = 1;
00206                 std::cout << "Now get the third column of the Jacobian forwardAD(1, x1)"  << std::endl;
00207                 funVals = osinstance->forwardAD(1, x1);
00208                 for( kjl = 0; kjl < 3; kjl++){
00209                         std::cout << "forward 1 " << funVals[ kjl] << std::endl;
00210                 }
00211                 
00220                 x1[0] = 1;
00221                 x1[1] = 0;
00222                 x1[2] = 1;
00223                 // recalculate the forward sweep with the new x1 vector
00224                 funVals = osinstance->forwardAD(1, x1);
00225                 std::vector<double> x2( n);
00226                 x2[0] = 0;
00227                 x2[1] = 0;
00228                 x2[2] = 0;
00229                 std::cout << "Now calcuate forwardAD(2, x2)"  << std::endl;
00230                 funVals = osinstance->forwardAD(2, x2);
00231                 for( kjl = 0; kjl < 3; kjl++){
00232                         std::cout << "forward 2 " << funVals[ kjl] << std::endl;
00233                 }                
00234                 
00242                 std::vector<double> vlambda(3);
00243                 vlambda[0] = 0;
00244                 vlambda[1] = 0;
00245                 vlambda[2] = 1;
00246                 // reverse sweep to get third row of Jacobian 
00247                 std::cout << "Now get the third row of the Jacobian reverseAD(1, vlambda)"  << std::endl;
00248                 osinstance->forwardAD(0, x0);
00249                 funVals = osinstance->reverseAD(1, vlambda);
00250                 for( kjl = 0; kjl < 3; kjl++){
00251                         std::cout << "reverse 1 " << funVals[ kjl] << std::endl;
00252                 }
00253                 // now get the Hessian of the Lagrangian of objective and 
00254                 // with the following multipliers
00255                 vlambda[0] = 1;
00256                 vlambda[1] = 2;
00257                 vlambda[2] = 1;
00263                 x1[0] = 0;
00264                 x1[1] = 0;
00265                 x1[2] = 1;
00266                 funVals = osinstance->forwardAD(1, x1);          
00267                  
00278                 dfunVals = osinstance->reverseAD(2, vlambda);
00279                 // get the first partials of the Lagrangian
00280                 std::cout << "Here are the first partials of the Lagrangain" << std::endl;
00281                 for(int kjl = 0; kjl <= 4; kjl+=2){
00282                         std::cout << dfunVals[ kjl] << std::endl;
00283                 }
00289                 std::cout << "Here is the third row (column) of Hessian of Lagrangian" << std::endl;
00290                 for(int kjl = 1; kjl <= 5; kjl+=2){
00291                         std::cout << dfunVals[ kjl] << std::endl;
00292                 }
00298                 double* x = new double[4]; //primal variables
00299                 double* z = new double[2]; //Lagrange multipliers on constraints
00300                 double* w = new double[1]; //Lagrange multiplier on objective
00301                 x[ 0] = 1;    // primal variable 0
00302                 x[ 1] = 5;    // primal variable 1
00303                 x[ 2] = 10;   // primal variable 2
00304                 x[ 3] = 5;    // primal variable 3
00305                 z[ 0] = 2;    // Lagrange multiplier on constraint 0
00306                 z[ 1] = 1;    // Lagrange multiplier on constraint 1
00307                 w[ 0] = 1;    // Lagrange multiplier on the objective function
00308                 
00338                 osinstance->bUseExpTreeForFunEval = false;
00339                 std::cout << "Calculate objective, idx = -1"  << std::endl;                     
00340                 std::cout << "obj value = " << osinstance->calculateFunctionValue(-1, x,  true) << std::endl;
00341                 
00342                 std::cout << "Calculate  first constraint, idx = 0"  << std::endl;                      
00343                 std::cout << "constraint index 0 value = " << osinstance->calculateFunctionValue(0, x,  true) << std::endl;
00344                  
00345                 std::cout << "Now use calculateAllConstraintFunctionValues"  << std::endl;                      
00346                 conVals = osinstance->calculateAllConstraintFunctionValues(x, true);
00347                 for( idx = 0; idx < osinstance->getConstraintNumber(); idx++){
00348                         std::cout << "CONSTRAINT FUNCTION INDEX = " <<  idx << "  CONSTRAINT FUNCTION VALUE =  "  << *(conVals + idx) << std::endl;
00349                 }
00350                 //
00351                 std::cout << "Now use calculateAllObjectiveFunctionValues"  << std::endl;       
00352                 objVals = osinstance->calculateAllObjectiveFunctionValues( x, NULL, NULL, true, 0);
00353                 for( idx = 0; idx < osinstance->getObjectiveNumber(); idx++){
00354                         std::cout << "OBJECTIVE FUNCTION  INDEX = " << idx <<  "  OBJECTIVE FUNCTION VALUE = "  << *(objVals + idx) << std::endl;
00355                 }
00356                 ok = CheckFunctionValues( conVals, *objVals, x[ 0], x[1], x[2], x[3],  z[0], z[1], w[0] );
00357                 if( ok == 0){
00358                         std::cout << "FAILED CHECKING FUNCTION VALUES TEST" << std::endl;
00359                         return 0;
00360                 }
00361                 else{
00362                         std::cout << "PASSED CHECKING FUNCTION VALUES TEST" << std::endl;
00363                 }
00364                 
00372                 SparseJacobianMatrix *sparseJac;
00373                 sparseJac = osinstance->getJacobianSparsityPattern();           
00374                 // print out just the sparsity pattern
00375                 std::cout << "JACOBIAN SPARSITY PATTERN"   << std::endl;
00376                 for(idx = 0; idx < sparseJac->startSize; idx++){
00377                         // some solvers (e.g. reduced gradient solvers) may want to know which values of the
00378                         // Jacobian matrix are constant, i.e. linear, sparseJac->conVals is the number of constant
00379                         // terms in the gradient for each rowt, the first conVals terms are constant, when getting
00380                         std::cout << "number constant terms in constraint "   <<  idx << " is " 
00381                         << *(sparseJac->conVals + idx)  << std::endl;
00382                         for(k = *(sparseJac->starts + idx); k < *(sparseJac->starts + idx + 1); k++){
00383                                 std::cout << "row idx = " << idx <<  "  col idx = "<< *(sparseJac->indexes + k) << std::endl;
00384                         }
00385                 }       
00386                 
00387                 SparseHessianMatrix *sparseHessian;
00388                 // the Hessian test
00389                 // get the sparsity pattern -- many solvers want to initialize with just the sparsity
00390                 std::cout << "GET LAGRANGIAN HESSIAN SPARSITY PATTERN"   << std::endl;
00391                 sparseHessian = osinstance->getLagrangianHessianSparsityPattern( );
00392                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){
00393                         std::cout <<  "Row Index = " << *(sparseHessian->hessRowIdx + idx) ;
00394                         std::cout <<  "  Column Index = " << *(sparseHessian->hessColIdx + idx) << std::endl;
00395                 }                
00396                 
00412                 double *objGrad;
00413                 std::cout << "OBJECTIVE FUNCTION GRADIENT"   << std::endl;
00414                 // in our implementation the objective function is a dense gradient
00415                 objGrad = osinstance->calculateObjectiveFunctionGradient( x, NULL, NULL,  -1, false, 1);
00416                 for(idx = 0; idx < osinstance->getVariableNumber(); idx++){
00417                         std::cout << "col idxx = " << idx << "  value =  " << *(objGrad + idx)  << std::endl;
00418                 }
00419                 std::cout << "CONSTRAINT JACOBIAN MATRIX"   << std::endl;
00420                 // now make the gradient calculations and fill in the sparse Jacobian matrix
00421                 sparseJac = osinstance->calculateAllConstraintFunctionGradients( x, NULL, NULL,  false, 1);
00422                 for(idx = 0; idx < osinstance->getConstraintNumber(); idx++){
00423                         for(k = *(sparseJac->starts + idx); k < *(sparseJac->starts + idx + 1); k++){
00424                                 std::cout << "row idx = " << idx <<  "  col idx = "<< *(sparseJac->indexes + k)
00425                                 << " value = " << *(sparseJac->values + k) << std::endl;
00426                         }
00427                 }
00428                 ok = CheckGradientValues( sparseJac, objGrad, x[ 0], x[1], x[2], x[3], z[0], z[1], w[0] );
00429                 if( ok == 0){
00430                         std::cout << "FAILED THE GRADIENT TEST" << std::endl;
00431                         return 0;
00432                 }
00433                 else{
00434                         std::cout << "PASSED THE GRADIENT TEST" << std::endl;
00435                 }                                         
00445                 //first iteration 
00446                 std::cout << "GET LAGRANGIAN HESSIAN FIRST TIME"   << std::endl;
00447                 sparseHessian = osinstance->calculateLagrangianHessian( x, w,  z,  false, 2);
00448                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){
00449                         std::cout << "row idx = " << *(sparseHessian->hessRowIdx + idx) <<  
00450                         "  col idx = "<< *(sparseHessian->hessColIdx + idx)
00451                         << " value = " << *(sparseHessian->hessValues + idx) << std::endl;
00452                 }
00453                 ok = CheckHessianUpper( sparseHessian, x[0],  x[1], x[2],  x[3],  z[0], z[1], w[0]);
00454                 if( ok == 0){
00455                         std::cout << "FAILED THE FIRST HESSIAN TEST" << std::endl;
00456                         return 0; 
00457                 }
00458                 else{
00459                         std::cout << "PASSED THE FIRST HESSIAN TEST" << std::endl;
00460                 }
00461                 // now change an x value, we don't rebuild the tree, however new_x 
00462                 // must be set to true
00463                 x[0] = 5;
00464                 
00465                 
00466                 
00467                 std::cout << "NOW GET LAGRANGIAN HESSIAN SECOND TIME FOR x[0] = 5"   << std::endl;
00468                 
00476                 osinstance->calculateAllObjectiveFunctionValues( x, w, z, true, 2);
00490                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){
00491                         std::cout << "row idx = " << *(sparseHessian->hessRowIdx + idx) <<  
00492                         "  col idx = "<< *(sparseHessian->hessColIdx + idx)
00493                         << " value = " << *(sparseHessian->hessValues + idx) << std::endl;
00494                 }
00495                 ok = CheckHessianUpper( sparseHessian , x[0],  x[1], x[2], x[3],  z[0], z[1], w[0] );
00496                 if( ok == 0){
00497                         std::cout << "FAILED THE SECOND HESSIAN TEST" << std::endl;
00498                         return 0;
00499                 }
00500                 else{
00501                         std::cout << "PASSED THE SECOND HESSIAN TEST" << std::endl  << std::endl ;
00502                 }
00506                 std::cout << "HERE IS ROW 1 OF JACOBIAN MATRIX" << std::endl; 
00507                 idx = 1;
00508                 for(k = *(sparseJac->starts + idx); k < *(sparseJac->starts + idx + 1); k++){
00509                         std::cout << "row idx = " << idx <<  "  col idx = "<< *(sparseJac->indexes + k)
00510                                 << " value = " << *(sparseJac->values + k) << std::endl;
00511                 }
00512                 std::cout << std::endl; 
00519                 sparseHessian = osinstance->calculateHessian(x, 1, true);
00520                 std::cout << "HERE IS ROW 1 HESSIAN MATRIX" << std::endl;
00521                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){
00522                         std::cout << "row idx = " << *(sparseHessian->hessRowIdx + idx) <<  
00523                         "  col idx = "<< *(sparseHessian->hessColIdx + idx)
00524                         << " value = " << *(sparseHessian->hessValues + idx) << std::endl;
00525                 }
00526                 //
00527                 // adjust the Lagrange multipliers to correspond to finding Hessian of constraint 1
00528                 z[ 0] = 0;  // Lagrange multiplier on constraint 0
00529                 z[ 1] = 1;  // Lagrange multiplier on constraint 1
00530                 w[ 0] = 0;  // Lagrange multiplier on the objective function
00531                 ok = CheckHessianUpper( sparseHessian , x[0],  x[1], x[2], x[3], z[0], z[1], w[0] );
00532                 if( ok == 0){
00533                         std::cout << "FAILED THE THIRD HESSIAN TEST" << std::endl;
00534                         return 0;
00535                 }
00536                 else{
00537                         std::cout << "PASSED THE THIRD HESSIAN TEST" << std::endl  << std::endl ;
00538                 }
00539                 //set x[0] back to its original value of 1
00540                 x[ 0] = 1;
00541                 //return 0;     
00542                 //
00543                 //
00544                 // now work directly with the CppAD package instead of OSInstance API
00545                 //
00546                 n = 4;
00547                 m = 3;
00548                 CppADvector< AD<double> > X(n);
00549                 CppADvector< AD<double> > Y(m);
00550                 X[0] = 5;
00551                 X[1] = 5;
00552                 X[2] = 0;
00553                 X[3] = 1;
00554                 // declare independent variables and start tape recording
00555                 std::cout << "Start Taping" << std::endl;
00556                 CppAD::Independent( X);
00557                 // range space vector 
00558                 // we include the constant terms in the CppAD functions
00559                 Y[ 0] =  CppAD::pow(X[0], 2) + 9*X[1];
00560                 Y[ 1] =  33 - 105 + 1.37*X[1] + 2*X[3] + 5*X[1] ;
00561                 Y[ 2] =  log(X[0]*X[3]) + 7*X[2] ;
00562                 // create f: x -> y and stop tape recording
00563                 CppAD::ADFun<double> f(X, Y); 
00564                 std::cout << "Stop Taping" << std::endl;
00565                 // get function values
00566                 std::vector<double> x_vec( n);
00567                 x_vec[ 0] = x[ 0];
00568                 x_vec[ 1] = x[ 1];
00569                 x_vec[ 2] = x[ 2];
00570                 x_vec[ 3] = x[ 3];
00571                 funVals = f.Forward(0, x_vec);
00572                 conVals[ 0] = funVals[ 1];
00573                 std::cout << "conVals[ 0] = " << conVals[ 0] << std::endl;
00574                 conVals[ 1] = funVals[ 2];
00575                 std::cout << "conVals[ 1] = " << conVals[ 1] << std::endl;
00576                 objVals[ 0] = funVals[ 0];
00577                 std::cout << "objVals[ 0] = " << objVals[ 0] << std::endl;
00578                 ok = CheckFunctionValues( conVals, funVals[ 0], x[ 0], x[1], x[2], x[3], z[0], z[1], w[0] );
00579                 if( ok == 0){
00580                         std::cout << "FAILED CHECKING FUNCTION VALUES TEST" << std::endl;
00581                         return 0;
00582                 }
00583                 else{
00584                         std::cout << "PASSED CHECKING FUNCTION VALUES TEST" << std::endl;
00585                 }
00586                 // now get gradient and Hessian
00587                 // first define and initialze unit vector vector
00588                 sparseJac = osinstance->getJacobianSparsityPattern();
00589                 std::vector<double> unit_col_vec( n);
00590                 std::vector<double> lagMultipliers( m); 
00591                 std::vector<double> gradVals( m);
00592                 lagMultipliers[ 0] = w[ 0];
00593                 lagMultipliers[ 1] = z[ 0];
00594                 lagMultipliers[ 2] = z[ 1];     
00595                 unsigned int index, kj;
00596                 return 0;
00597                 for(index = 0; index < n; index++){
00598                         unit_col_vec[ index] = 0;
00599                 }       
00600                 for(index = 0; index < n; index++){
00601                         unit_col_vec[ index] = 1;
00602                         // calculate column i of the Jacobian matrix
00603                         gradVals = f.Forward(1, unit_col_vec);
00604                         unit_col_vec[ index] = 0;
00605                         // get the nonzero gradient values in constraint k
00606                         for(kj = 0; kj < m; kj++){
00607                                 std::cout << "variable " << index << "  row " << kj << "  gradient value" << std::endl;
00608                                 std::cout << "gradient value = " << gradVals[ kj] << std::endl; 
00609                         }
00610                         // get row i of the Lagrangian function!!!
00611                         f.Reverse(2, lagMultipliers);
00612                 }
00613                 // done with CppAD test 
00614                 // do garbage collection
00615                 delete osilreader;
00616                 osilreader = NULL;
00617                 std::cout << "OSILREADER DELETED" << std::endl; 
00618                 delete[] conVals;
00619                 delete[] objVals;               
00620         }
00621         catch(const ErrorClass& eclass){
00622                 std::cout << eclass.errormsg << std::endl;
00623         } 
00624 
00625         //
00626         {
00627                 //checking CppAD power, another illustration of CppAD
00628                 size_t n  = 2;
00629                 double x0 = 4;
00630                 double x1 = .5;
00631                 CppADvector< AD<double> > x(n);
00632             x[0]      = x0;
00633             x[1]      = x1;
00634              // declare independent variables and start tape recording
00635              CppAD::Independent(x);
00636              // range space vector 
00637              size_t m = 1;
00638              CppADvector< AD<double> > y(m);
00639              y[0] = CppAD::pow(x[0], x[1]);
00640              // create f: x -> y and stop tape recording
00641              CppAD::ADFun<double> f(x, y); 
00642              // check value 
00643              double check = std::pow(x0, x1);
00644              check = CppAD::pow(x0, x1);
00645              std::cout << "check " <<  y[ 0] << std::endl;
00646              ok &= NearEqual(y[0] , check,  1e-10 , 1e-10);
00647              // forward computation of first partial w.r.t. x[0]
00648              std::vector<double> dx(n);
00649              std::vector<double> dy(m);
00650              dx[0] = 4.;
00651              dx[1] = 1/2.;
00652              dy    = f.Forward(1, dx);
00653              std::cout << "dy =  " <<  dy[ 0] << std::endl;
00654              check = x1 * std::pow(x0, x1-1.);
00655              ok   &= NearEqual(dy[0], check, 1e-10, 1e-10);
00656         }
00657         return 0;
00658 }// end main program
00659 
00660 bool CheckFunctionValues( double *conVals, double objValue,
00661         double x0, double x1, double x2, double x3, double z0, double z1, double w ){
00662         using CppAD::NearEqual;
00663         bool ok  = true;
00664         double checkObj = x0*x0 + 9*x1;
00665         std::cout  << "checkObj = " << checkObj << std::endl;
00666         ok &= NearEqual(objValue, checkObj, 1e-10, 1e-10); 
00667         double checkCon0 = 33. - 105. + 1.37*x1 + 2*x3 + 5*x1;
00668         std::cout  << "checkCon0 = " << checkCon0 << std::endl;
00669         ok &= NearEqual(*(conVals + 0), checkCon0, 1e-10, 1e-10);
00670         double checkCon1 = log(x0*x3) + 7*x2;
00671         std::cout  << "checkCon1 = " << checkCon1 << std::endl;
00672         ok &= NearEqual( *(conVals + 1), checkCon1, 1e-10, 1e-10);
00673         return ok;
00674 }//CheckFunctionValues
00675 //
00676 //
00677 bool CheckGradientValues( SparseJacobianMatrix *sparseJac, double *objGrad,
00678         double x0, double x1, double x2, double x3, double y0, double y1, double w ){
00679         using CppAD::NearEqual;
00680         bool ok  = true;
00681         // first the objective function gradient
00682         double checkObjPartial0 = 2*x0;
00683         ok &= NearEqual( *(objGrad + 0), checkObjPartial0, 1e-10, 1e-10); 
00684         double checkObjPartial1 = 9;
00685         ok &= NearEqual( *(objGrad + 1), checkObjPartial1, 1e-10, 1e-10);       
00686         double checkObjPartial2 = 0;
00687         ok &= NearEqual( *(objGrad + 2), checkObjPartial2, 1e-10, 1e-10); 
00688         // get the constrating gradient
00689         // row 0 gradient -- there are nonzero partials for variables 1 and 2
00690         double checkCon0Partial1 = 1.37 + 5.0;
00691         ok &= NearEqual( *(sparseJac->values + 0), checkCon0Partial1, 1e-10, 1e-10);    
00692         double checkCon0Partial3 = 2.;
00693         ok &= NearEqual( *(sparseJac->values + 1), checkCon0Partial3, 1e-10, 1e-10); 
00694         // row 1 gradient -- there are nonzero partials for variables 0 and 2
00695         double checkCon1Partial2 = 7;
00696         ok &= NearEqual( *(sparseJac->values + 2), checkCon1Partial2, 1e-10, 1e-10);    
00697         double checkCon1Partial0 = 1./x0;
00698         ok &= NearEqual( *(sparseJac->values + 3), checkCon1Partial0, 1e-10, 1e-10);    
00699         double checkCon1Partial3 = 1./x3;
00700         ok &= NearEqual( *(sparseJac->values + 4), checkCon1Partial3, 1e-10, 1e-10); 
00701         return ok;
00702 }//CheckGradientValues
00703 //
00704 bool CheckHessianUpper( SparseHessianMatrix *sparseHessian , 
00705         double x0, double x1, double x2, double x3, double z0, double z1, double w ){
00706         using CppAD::NearEqual;
00707         bool ok  = true;
00708         int hessValuesIdx = 0;
00709         //assert( sparseHessian->hessDimension = n * (n + 1) /2)
00710         // L_00 = 2 * w - z1 / ( x0 * x0 )
00711         double check = 2. * w - z1 / (x0 * x0);
00712         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), check, 1e-10, 1e-10); 
00713         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), 0., 1e-10, 1e-10);
00714         if(ok == false) std::cout << "FAILED TWO" << std::endl;
00715         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), 0., 1e-10, 1e-10);
00716         if(ok == false) std::cout << "FAILED THREE" << std::endl;
00717         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), 0., 1e-10, 1e-10);
00718         if(ok == false) std::cout << "FAILED FOUR" << std::endl;
00719         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), 0., 1e-10, 1e-10);
00720         if(ok == false) std::cout << "FAILED FIVE" << std::endl;
00721         // L_22 = - z1 / (x3 * x3)
00722         check = - z1 / (x3 * x3);
00723         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), check, 1e-10, 1e-1);
00724         if(ok == false) std::cout << "FAILED SIX" << std::endl;
00725         return ok;
00726 }//CheckHessianUpper

Generated on Thu May 15 22:15:04 2008 by  doxygen 1.4.7