- Main Page
- Namespaces
- Classes
- Files
- Related Pages

- File List
- File Members

# /home/coin/SVN-release/OS-2.3.5/ApplicationTemplates/algorithmicDiff/OSAlgorithmicDiffTest.cpp

Go to the documentation of this file.

    00001 /* $Id: OSAlgorithmicDiffTest.cpp 2698 2009-06-09 04:14:07Z kmartin $ */
    00075 #include <cstddef>
    00076 #include <cstdlib>
    00077 #include <cctype>
    00078 #include <cassert>
    00079 #include <stack>
    00080 #include <cppad/cppad.hpp>
    00081 #include <iostream>
    00082
    00083 // CoinHelperFunctions <has cmath>
    00084 #include "CoinHelperFunctions.hpp"
    00085 #include "OSInstance.h"
    00086 #include "OSiLWriter.h"
    00087 #include "OSParameters.h"
    00088 #include "OSnLNode.h"
    00089 #include "OSErrorClass.h"
    00090 #include "OSFileUtil.h"
    00091 #include "OSiLReader.h"
    00092 #include "OSInstance.h"
    00093 #include "OSExpressionTree.h"
    00094 #include "OSnLNode.h"
    00095 #include "OSDataStructures.h"
    00096
    00097
    00098
    00099 #include <vector>
    00100 #include <map>
    00101 #include<string>
    00102
    00103
    00104 int  main(){
    00105         WindowsErrorPopupBlocker();
    00106         using std::cout;
    00107         using std::endl;
    00108         using CppAD::AD;
    00109         //using CppAD::NearEqual;
    00110         using CppAD::vector;
    00111         std::cout.precision(12);
    00112         // error checking functions
    00113         bool CheckFunctionValues( double *conVals, double objValue,
    00114                 double x0, double x1, double x2, double x3, double y0, double y1, double z );
    00115         bool CheckHessianUpper( SparseHessianMatrix *sparseHessian,
    00116                 double x0, double x1, double x2, double x3, double y0, double y1, double z );
    00117         bool CheckGradientValues( SparseJacobianMatrix *sparseJac, double *objGrad,
    00118                 double x0, double x1, double x2, double x3, double y0, double y1, double z );
    00119         bool ok = true;
    00120         int k, idx;
    00121         //
    00122         // get the problem data
    00123         //
    00124         FileUtil *fileUtil = NULL;
    00125         std::string osilFileName;
    00126         std::string osil;
    00127         // get the input file
    00128         const char dirsep =  CoinFindDirSeparator();
    00129         // Set directory containing mps data files.
    00130         std::string dataDir;
    00131     dataDir = dirsep == '/' ? "../data/" : "..\\data\\";
    00132         //osilFileName =  dataDir + "HS071_NLP.osil";
    00133         osilFileName =  dataDir   + "CppADTestLag.osil";
    00134         std::cout  << "osilFileName  =  " << osilFileName << std::endl;
    00135         fileUtil = new FileUtil();
    00136         osil = fileUtil->getFileAsString( &osilFileName[0]);
    00137         //
    00138         // create OSReader and OSInstance objects
    00139         OSiLReader *osilreader = NULL;
    00140         OSInstance *osinstance = NULL;
    00141         // create reader, parse the OSiL, and generate the OSInstance object
    00142         try{
    00143                 // a counter
    00144                 int kjl;
    00145                 osilreader = new OSiLReader();
    00146                 osinstance = osilreader->readOSiL( osil);
    00147                 std::vector<double> funVals(3);
    00148                 std::vector<double> dfunVals(6);
    00149                 double *conVals = NULL;
    00150                 //conVals = new double[ 2];
    00151                 double *objVals = NULL;
    00152                 //objVals = new double[ 1];
    00153
    00167                 //
    00168                 // first initialize the nonlinear structures for call backs
    00169                 std::cout << "Initialize Nonlinear Structures" << std::endl;
    00170                 osinstance->initForAlgDiff( );
    00171                 //osinstance->getJacobianSparsityPattern( );
    00172
    00173
    00174
    00180                 std::map<int, int> varIndexMap;
    00181                 std::map<int, int>::iterator posVarIndexMap;
    00182                 varIndexMap = osinstance->getAllNonlinearVariablesIndexMap( );
    00187                 for(posVarIndexMap = varIndexMap.begin(); posVarIndexMap != varIndexMap.end(); ++posVarIndexMap){
    00188                                 std::cout <<  "Variable Index = "   << posVarIndexMap->first  << std::endl ;
    00189                 }
    00190                 std::cout << "Number of nonlinear variables =  " << varIndexMap.size() << std::endl;
    00191
    00192                 //
    00193
    00194                 // get the number of nonlinear terms
    00195
    00196                 int mm = osinstance->getNumberOfNonlinearExpressionTreeModIndexes();
    00197
    00198                 int jj;
    00199
    00200                 for(jj = 0; jj < mm; jj++){
    00201                         std::cout << osinstance->getNonlinearExpressionTreeModIndexes()[ jj] << std::endl;
    00202                 }
    00203
    00204                 std::cout << "Number of unique nonlinear terms =  " <<  mm << std::endl;
    00205                 //return 0;
    00206
    00207
    00208
    00209                 // domain space vector
    00210                 size_t n  = varIndexMap.size(); // three variables
    00211                 // range space vector
    00212                 size_t m = 3; // Lagrangian has an objective and two constraints
    00213
    00214                 std::vector<double> x0( n);
    00219                 x0[0] = 1; // the value for variable x0
    00220                 x0[1] = 5; // the value for variable x1
    00221                 x0[2] = 5; // the value for variable x3
    00222                 std::cout << "CALL forward" << std::endl;
    00223                 funVals = osinstance->forwardAD(0, x0);
    00224                 for( kjl = 0; kjl < 3; kjl++){
    00225                         std::cout << "forward 0 " << funVals[ kjl] << std::endl;
    00226                 }
    00227                 // get the third column of the Jacobian from a forward sweep
    00228                 std::vector<double> x1( n);
    00229                 x1[0] = 0;
    00230                 x1[1] = 0;
    00231                 x1[2] = 1;
    00232                 std::cout << "Now get the third column of the Jacobian forwardAD(1, x1)"  << std::endl;
    00233                 funVals = osinstance->forwardAD(1, x1);
    00234                 for( kjl = 0; kjl < 3; kjl++){
    00235                         std::cout << "forward 1 " << funVals[ kjl] << std::endl;
    00236                 }
    00237
    00246                 x1[0] = 1;
    00247                 x1[1] = 0;
    00248                 x1[2] = 1;
    00249                 // recalculate the forward sweep with the new x1 vector
    00250                 funVals = osinstance->forwardAD(1, x1);
    00251                 std::vector<double> x2( n);
    00252                 x2[0] = 0;
    00253                 x2[1] = 0;
    00254                 x2[2] = 0;
    00255                 std::cout << "Now calcuate forwardAD(2, x2)"  << std::endl;
    00256                 funVals = osinstance->forwardAD(2, x2);
    00257                 for( kjl = 0; kjl < 3; kjl++){
    00258                         std::cout << "forward 2 " << funVals[ kjl] << std::endl;
    00259                 }
    00260
    00268                 std::vector<double> vlambda(3);
    00269                 vlambda[0] = 0;
    00270                 vlambda[1] = 0;
    00271                 vlambda[2] = 1;
    00272                 // reverse sweep to get third row of Jacobian
    00273                 std::cout << "Now get the third row of the Jacobian reverseAD(1, vlambda)"  << std::endl;
    00274                 osinstance->forwardAD(0, x0);
    00275                 funVals = osinstance->reverseAD(1, vlambda);
    00276                 for( kjl = 0; kjl < 3; kjl++){
    00277                         std::cout << "reverse 1 " << funVals[ kjl] << std::endl;
    00278                 }
    00279                 // now get the Hessian of the Lagrangian of objective and
    00280                 // with the following multipliers
    00281                 vlambda[0] = 1;
    00282                 vlambda[1] = 2;
    00283                 vlambda[2] = 1;
    00289                 x1[0] = 0;
    00290                 x1[1] = 0;
    00291                 x1[2] = 1;
    00292                 funVals = osinstance->forwardAD(1, x1);
    00293
    00304                 dfunVals = osinstance->reverseAD(2, vlambda);
    00305                 // get the first partials of the Lagrangian
    00306                 std::cout << "Here are the first partials of the Lagrangain" << std::endl;
    00307                 for(int kjl = 0; kjl <= 4; kjl+=2){
    00308                         std::cout << dfunVals[ kjl] << std::endl;
    00309                 }
    00315                 std::cout << "Here is the third row (column) of Hessian of Lagrangian" << std::endl;
    00316                 for(int kjl = 1; kjl <= 5; kjl+=2){
    00317                         std::cout << dfunVals[ kjl] << std::endl;
    00318                 }
    00324                 double* x = new double[4]; //primal variables
    00325                 double* z = new double[2]; //Lagrange multipliers on constraints
    00326                 double* w = new double[1]; //Lagrange multiplier on objective
    00327                 x[ 0] = 1;    // primal variable 0
    00328                 x[ 1] = 5;    // primal variable 1
    00329                 x[ 2] = 10;   // primal variable 2
    00330                 x[ 3] = 5;    // primal variable 3
    00331                 z[ 0] = 2;    // Lagrange multiplier on constraint 0
    00332                 z[ 1] = 1;    // Lagrange multiplier on constraint 1
    00333                 w[ 0] = 1;    // Lagrange multiplier on the objective function
    00334
    00364                 osinstance->bUseExpTreeForFunEval = false;
    00365                 std::cout << "Calculate objective, idx = -1"  << std::endl;
    00366                 std::cout << "obj value = " << osinstance->calculateFunctionValue(-1, x,  true) << std::endl;
    00367
    00368                 std::cout << "Calculate  first constraint, idx = 0"  << std::endl;
    00369                 std::cout << "constraint index 0 value = " << osinstance->calculateFunctionValue(0, x,  true) << std::endl;
    00370
    00371                 std::cout << "Now use calculateAllConstraintFunctionValues"  << std::endl;
    00372                 conVals = osinstance->calculateAllConstraintFunctionValues(x, true);
    00373                 for( idx = 0; idx < osinstance->getConstraintNumber(); idx++){
    00374                         std::cout << "CONSTRAINT FUNCTION INDEX = " <<  idx << "  CONSTRAINT FUNCTION VALUE =  "  << *(conVals + idx) << std::endl;
    00375                 }
    00376                 //
    00377                 std::cout << "Now use calculateAllObjectiveFunctionValues"  << std::endl;
    00378                 objVals = osinstance->calculateAllObjectiveFunctionValues( x, NULL, NULL, true, 0);
    00379                 for( idx = 0; idx < osinstance->getObjectiveNumber(); idx++){
    00380                         std::cout << "OBJECTIVE FUNCTION  INDEX = " << idx <<  "  OBJECTIVE FUNCTION VALUE = "  << *(objVals + idx) << std::endl;
    00381                 }
    00382                 ok = CheckFunctionValues( conVals, *objVals, x[ 0], x[1], x[2], x[3],  z[0], z[1], w[0] );
    00383                 if( ok == 0){
    00384                         std::cout << "FAILED CHECKING FUNCTION VALUES TEST" << std::endl;
    00385                         return 0;
    00386                 }
    00387                 else{
    00388                         std::cout << "PASSED CHECKING FUNCTION VALUES TEST" << std::endl;
    00389                 }
    00390
    00397                 std::cout << "PLACE CALL TO JACOBIAN SPARSITY PATTERN"   << std::endl;
    00398                 SparseJacobianMatrix *sparseJac;
    00399                 sparseJac = osinstance->getJacobianSparsityPattern();
    00400                 // print out just the sparsity pattern
    00401                 std::cout << "JACOBIAN SPARSITY PATTERN"   << std::endl;
    00402                 std::cout << "JACOBIAN START SIZE "   <<  sparseJac->startSize << std::endl;
    00403                 for(idx = 0; idx < osinstance->getConstraintNumber(); idx++){
    00404                         // some solvers (e.g. reduced gradient solvers) may want to know which values of the
    00405                         // Jacobian matrix are constant, i.e. linear, sparseJac->conVals is the number of constant
    00406                         // terms in the gradient for each rowt, the first conVals terms are constant, when getting
    00407                         std::cout << "number constant terms in constraint "   <<  idx << " is "
    00408                         << *(sparseJac->conVals + idx)  << std::endl;
    00409                         for(k = *(sparseJac->starts + idx); k < *(sparseJac->starts + idx + 1); k++){
    00410                                 std::cout << "row idx = " << idx <<  "  col idx = "<< *(sparseJac->indexes + k) << std::endl;
    00411                         }
    00412                 }
    00413
    00414                 SparseHessianMatrix *sparseHessian;
    00415                 // the Hessian test
    00416                 // get the sparsity pattern -- many solvers want to initialize with just the sparsity
    00417                 std::cout << "GET LAGRANGIAN HESSIAN SPARSITY PATTERN"   << std::endl;
    00418                 sparseHessian = osinstance->getLagrangianHessianSparsityPattern( );
    00419                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){
    00420                         std::cout <<  "Row Index = " << *(sparseHessian->hessRowIdx + idx) ;
    00421                         std::cout <<  "  Column Index = " << *(sparseHessian->hessColIdx + idx) << std::endl;
    00422                 }
    00423
    00439                 double *objGrad;
    00440                 std::cout << "OBJECTIVE FUNCTION GRADIENT"   << std::endl;
    00441                 // in our implementation the objective function is a dense gradient
    00442                 objGrad = osinstance->calculateObjectiveFunctionGradient( x, NULL, NULL,  -1, false, 1);
    00443                 for(idx = 0; idx < osinstance->getVariableNumber(); idx++){
    00444                         std::cout << "col idxx = " << idx << "  value =  " << *(objGrad + idx)  << std::endl;
    00445                 }
    00446                 std::cout << "CONSTRAINT JACOBIAN MATRIX"   << std::endl;
    00447                 // now make the gradient calculations and fill in the sparse Jacobian matrix
    00448                 sparseJac = osinstance->calculateAllConstraintFunctionGradients( x, NULL, NULL,  false, 1);
    00449                 for(idx = 0; idx < osinstance->getConstraintNumber(); idx++){
    00450                         for(k = *(sparseJac->starts + idx); k < *(sparseJac->starts + idx + 1); k++){
    00451                                 std::cout << "row idx = " << idx <<  "  col idx = "<< *(sparseJac->indexes + k)
    00452                                 << " value = " << *(sparseJac->values + k) << std::endl;
    00453                         }
    00454                 }
    00455                 ok = CheckGradientValues( sparseJac, objGrad, x[ 0], x[1], x[2], x[3], z[0], z[1], w[0] );
    00456                 if( ok == 0){
    00457                         std::cout << "FAILED THE GRADIENT TEST" << std::endl;
    00458                         return 0;
    00459                 }
    00460                 else{
    00461                         std::cout << "PASSED THE GRADIENT TEST" << std::endl;
    00462                 }
    00472                 //first iteration
    00473                 std::cout << "GET LAGRANGIAN HESSIAN FIRST TIME"   << std::endl;
    00474                 sparseHessian = osinstance->calculateLagrangianHessian( x, w,  z,  false, 2);
    00475                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){
    00476                         std::cout << "row idx = " << *(sparseHessian->hessRowIdx + idx) <<
    00477                         "  col idx = "<< *(sparseHessian->hessColIdx + idx)
    00478                         << " value = " << *(sparseHessian->hessValues + idx) << std::endl;
    00479                 }
    00480                 ok = CheckHessianUpper( sparseHessian, x[0],  x[1], x[2],  x[3],  z[0], z[1], w[0]);
    00481                 if( ok == 0){
    00482                         std::cout << "FAILED THE FIRST HESSIAN TEST" << std::endl;
    00483                         return 0;
    00484                 }
    00485                 else{
    00486                         std::cout << "PASSED THE FIRST HESSIAN TEST" << std::endl;
    00487                 }
    00488                 // now change an x value, we don't rebuild the tree, however new_x
    00489                 // must be set to true
    00490                 x[0] = 5;
    00491
    00492
    00493
    00494                 std::cout << "NOW GET LAGRANGIAN HESSIAN SECOND TIME FOR x[0] = 5"   << std::endl;
    00495
    00503                 osinstance->calculateAllObjectiveFunctionValues( x, w, z, true, 2);
    00517                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){
    00518                         std::cout << "row idx = " << *(sparseHessian->hessRowIdx + idx) <<
    00519                         "  col idx = "<< *(sparseHessian->hessColIdx + idx)
    00520                         << " value = " << *(sparseHessian->hessValues + idx) << std::endl;
    00521                 }
    00522                 ok = CheckHessianUpper( sparseHessian , x[0],  x[1], x[2], x[3],  z[0], z[1], w[0] );
    00523                 if( ok == 0){
    00524                         std::cout << "FAILED THE SECOND HESSIAN TEST" << std::endl;
    00525                         return 0;
    00526                 }
    00527                 else{
    00528                         std::cout << "PASSED THE SECOND HESSIAN TEST" << std::endl  << std::endl ;
    00529                 }
    00533                 std::cout << "HERE IS ROW 1 OF JACOBIAN MATRIX" << std::endl;
    00534                 idx = 1;
    00535                 for(k = *(sparseJac->starts + idx); k < *(sparseJac->starts + idx + 1); k++){
    00536                         std::cout << "row idx = " << idx <<  "  col idx = "<< *(sparseJac->indexes + k)
    00537                                 << " value = " << *(sparseJac->values + k) << std::endl;
    00538                 }
    00539                 std::cout << std::endl;
    00546                 sparseHessian = osinstance->calculateHessian(x, 1, true);
    00547                 std::cout << "HERE IS ROW 1 HESSIAN MATRIX" << std::endl;
    00548                 for(idx = 0; idx < sparseHessian->hessDimension; idx++){
    00549                         std::cout << "row idx = " << *(sparseHessian->hessRowIdx + idx) <<
    00550                         "  col idx = "<< *(sparseHessian->hessColIdx + idx)
    00551                         << " value = " << *(sparseHessian->hessValues + idx) << std::endl;
    00552                 }
    00553                 //
    00554                 // adjust the Lagrange multipliers to correspond to finding Hessian of constraint 1
    00555                 z[ 0] = 0;  // Lagrange multiplier on constraint 0
    00556                 z[ 1] = 1;  // Lagrange multiplier on constraint 1
    00557                 w[ 0] = 0;  // Lagrange multiplier on the objective function
    00558                 ok = CheckHessianUpper( sparseHessian , x[0],  x[1], x[2], x[3], z[0], z[1], w[0] );
    00559                 if( ok == 0){
    00560                         std::cout << "FAILED THE THIRD HESSIAN TEST" << std::endl;
    00561                         return 0;
    00562                 }
    00563                 else{
    00564                         std::cout << "PASSED THE THIRD HESSIAN TEST" << std::endl  << std::endl ;
    00565                 }
    00566                 //set x[0] back to its original value of 1
    00567                 x[ 0] = 1;
    00568                 //return 0;
    00569                 //
    00570                 //
    00571                 // now work directly with the CppAD package instead of OSInstance API
    00572                 //
    00573                 n = 4;
    00574                 m = 3;
    00575                 CppADvector< AD<double> > X(n);
    00576                 CppADvector< AD<double> > Y(m);
    00577                 X[0] = 5;
    00578                 X[1] = 5;
    00579                 X[2] = 0;
    00580                 X[3] = 1;
    00581                 // declare independent variables and start tape recording
    00582                 std::cout << "Start Taping" << std::endl;
    00583                 CppAD::Independent( X);
    00584                 // range space vector
    00585                 // we include the constant terms in the CppAD functions
    00586                 Y[ 0] =  CppAD::pow(X[0], 2) + 9*X[1];
    00587                 Y[ 1] =  33 - 105 + 1.37*X[1] + 2*X[3] + 5*X[1] ;
    00588                 Y[ 2] =  log(X[0]*X[3]) + 7*X[2] ;
    00589                 // create f: x -> y and stop tape recording
    00590                 CppAD::ADFun<double> f(X, Y);
    00591                 std::cout << "Stop Taping" << std::endl;
    00592                 // get function values
    00593                 std::vector<double> x_vec( n);
    00594                 x_vec[ 0] = x[ 0];
    00595                 x_vec[ 1] = x[ 1];
    00596                 x_vec[ 2] = x[ 2];
    00597                 x_vec[ 3] = x[ 3];
    00598                 funVals = f.Forward(0, x_vec);
    00599                 conVals[ 0] = funVals[ 1];
    00600                 std::cout << "conVals[ 0] = " << conVals[ 0] << std::endl;
    00601                 conVals[ 1] = funVals[ 2];
    00602                 std::cout << "conVals[ 1] = " << conVals[ 1] << std::endl;
    00603                 objVals[ 0] = funVals[ 0];
    00604                 std::cout << "objVals[ 0] = " << objVals[ 0] << std::endl;
    00605                 ok = CheckFunctionValues( conVals, funVals[ 0], x[ 0], x[1], x[2], x[3], z[0], z[1], w[0] );
    00606                 if( ok == 0){
    00607                         std::cout << "FAILED CHECKING FUNCTION VALUES TEST" << std::endl;
    00608                         return 0;
    00609                 }
    00610                 else{
    00611                         std::cout << "PASSED CHECKING FUNCTION VALUES TEST" << std::endl;
    00612                 }
    00613                 // now get gradient and Hessian
    00614                 // first define and initialze unit vector vector
    00615                 sparseJac = osinstance->getJacobianSparsityPattern();
    00616                 std::vector<double> unit_col_vec( n);
    00617                 std::vector<double> lagMultipliers( m);
    00618                 std::vector<double> gradVals( m);
    00619                 lagMultipliers[ 0] = w[ 0];
    00620                 lagMultipliers[ 1] = z[ 0];
    00621                 lagMultipliers[ 2] = z[ 1];
    00622                 unsigned int index, kj;
    00623                 //return 0;
    00624                 for(index = 0; index < n; index++){
    00625                         unit_col_vec[ index] = 0;
    00626                 }
    00627                 for(index = 0; index < n; index++){
    00628                         unit_col_vec[ index] = 1;
    00629                         // calculate column i of the Jacobian matrix
    00630                         gradVals = f.Forward(1, unit_col_vec);
    00631                         unit_col_vec[ index] = 0;
    00632                         // get the nonzero gradient values in constraint k
    00633                         for(kj = 0; kj < m; kj++){
    00634                                 std::cout << "variable " << index << "  row " << kj << "  gradient value" << std::endl;
    00635                                 std::cout << "gradient value = " << gradVals[ kj] << std::endl;
    00636                         }
    00637                         // get row i of the Lagrangian function!!!
    00638                         std::cout << "CALL f.Reverse -------" << std::endl;
    00639                         f.Reverse(2, lagMultipliers);
    00640                         std::cout << "FINISH CALL f.Reverse -------" << std::endl;
    00641                 }
    00642                 // done with CppAD test
    00643                 // do garbage collection
    00644                 delete osilreader;
    00645                 osilreader = NULL;
    00646                 std::cout << "OSILREADER DELETED" << std::endl;
    00647                 //delete[] conVals;
    00648                 //delete[] objVals;
    00649                 delete[] x;
    00650                 delete[] z;
    00651                 delete[] w;
    00652         }
    00653         catch(const ErrorClass& eclass){
    00654                 std::cout << eclass.errormsg << std::endl;
    00655         }
    00656
    00657         //
    00658         {
    00659                 //checking CppAD power, another illustration of CppAD
    00660                 size_t n  = 2;
    00661                 double x0 = 4;
    00662                 double x1 = .5;
    00663                 CppADvector< AD<double> > x(n);
    00664             x[0]      = x0;
    00665             x[1]      = x1;
    00666              // declare independent variables and start tape recording
    00667              CppAD::Independent(x);
    00668              // range space vector
    00669              size_t m = 1;
    00670              CppADvector< AD<double> > y(m);
    00671              y[0] = std::pow(x0, x1);
    00672              // create f: x -> y and stop tape recording
    00673              CppAD::ADFun<double> f(x, y);
    00674              // check value
    00675              double check = std::pow(x0, x1);
    00676              // forward computation of first partial w.r.t. x[0]
    00677              std::vector<double> dx(n);
    00678              std::vector<double> dy(m);
    00679              dx[0] = 4.;
    00680              dx[1] = 1/2.;
    00681              dy    = f.Forward(1, dx);
    00682              std::cout << "dy =  " <<  dy[ 0] << std::endl;
    00683              check = x1 * std::pow(x0, x1-1.);
    00684              //ok   &= NearEqual(dy[0], check, 1e-10, 1e-10);
    00685              ok = ( fabs(check - dy[0])/(fabs( check) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00686         }
    00687
    00688         {
    00689
    00690                 //checking CppAD sparsity features
    00691                 // domain space vector
    00692                 size_t n = 3;
    00693                 CPPAD_TEST_VECTOR< AD<double> > X(n);
    00694                 X[0] = 0.;
    00695                 X[1] = 1.;
    00696                 X[2] = 2.;
    00697                 // declare independent variables and start recording
    00698                 CppAD::Independent(X);
    00699                 // range space vector
    00700                 size_t m = 2;
    00701                 CPPAD_TEST_VECTOR< AD<double> > Y(m);
    00702                 Y[0] = CppAD::pow(X[0], 2)  + CppAD::pow(X[2], 2);
    00703                 Y[1] = -CppAD::pow(X[0], 2) + CppAD::pow(X[1], 2);
    00704                 // create f: X -> Y and stop tape recording
    00705                 CppAD::ADFun<double> f(X, Y);
    00706
    00707                 // sparsity pattern for the identity matrix
    00708                 std::vector<bool> r(n * n);
    00709                 size_t i, j;
    00710                 for(i = 0; i < n; i++) {
    00711                         for(j = 0; j < n; j++)
    00712                                 r[ i * n + j ] = false;
    00713                                 r[ i * n + i ] = true;
    00714                 }
    00715                 // compute sparsity pattern for J(x) = F^{(1)} (x)
    00716                 f.ForSparseJac(n, r);
    00720                 //std::vector<bool> s(m * m);
    00721                 //for(i = 0; i < m; i++){
    00722                 //      for(j = 0; j < m; j++)
    00723                 //              s[ i * m + j ] = false;
    00724                 //      s[ i * m + i ] = true;
    00725             // }
    00726             // sparsity pattern for F'(x)
    00727             // f.RevSparseJac(m, s);
    00731                 // compute sparsity pattern for H(x) = F_0^{(2)} (x)
    00732                 std::vector<bool> e( m);
    00733                 //Vector s(m);
    00734                 for(i = 0; i < m; i++)
    00735                 e[i] = false;
    00736                 e[ 0] = true;
    00737                 e[ 1] = false;
    00738                 std::vector<bool> h( n*n);
    00739                 //Vector h(n * n);
    00740                 std::cout << "Computing Sparse Hessian" << std::endl;
    00741                 h = f.RevSparseHes(n, e);
    00742                 for(i = 0; i < n; i++){
    00743                         std::cout << "Row " << i << "  of Hessian " << std::endl;
    00744                         for(j = 0; j < n; j++){
    00745                                 std::cout << h[ i*n + j] <<  "  ";
    00746                         }
    00747                         std::cout << std::endl;
    00748                 }
    00749         }
    00750         delete fileUtil;
    00751         return 0;
    00752 }// end main program
    00753
    00754 bool CheckFunctionValues( double *conVals, double objValue,
    00755         double x0, double x1, double x2, double x3, double z0, double z1, double w ){
    00756         using CppAD::NearEqual;
    00757         bool ok  = true;
    00758         double checkObj = x0*x0 + 9*x1;
    00759         std::cout  << "checkObj = " << checkObj << std::endl;
    00760         std::cout  << "objValue = " << objValue << std::endl;
    00761         //ok &= NearEqual(objValue, checkObj, 1e-10, 1e-10);
    00762         ok = ( fabs(checkObj - objValue )/(fabs( checkObj) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00763         double checkCon0 = 33. - 105. + 1.37*x1 + 2*x3 + 5*x1;
    00764         std::cout  << "checkCon0 = " << checkCon0 << std::endl;
    00765         std::cout  << "conVals = " << *(conVals + 0) << std::endl;
    00766         //ok &= NearEqual(*(conVals + 0), checkCon0, 1e-10, 1e-10);
    00767         ok = ( fabs(checkCon0 - *(conVals + 0) )/(fabs( checkCon0) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00768         double checkCon1 = log(x0*x3) + 7*x2;
    00769         std::cout  << "checkCon1 = " << checkCon1 << std::endl;
    00770         std::cout  << "conVals = " << *(conVals + 1) << std::endl;
    00771         //ok &= NearEqual( *(conVals + 1), checkCon1, 1e-10, 1e-10);
    00772         ok = ( fabs(checkCon1 - *(conVals + 1) )/(fabs( checkCon1) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00773         return ok;
    00774 }//CheckFunctionValues
    00775 //
    00776 //
    00777 bool CheckGradientValues( SparseJacobianMatrix *sparseJac, double *objGrad,
    00778         double x0, double x1, double x2, double x3, double y0, double y1, double w ){
    00779         using CppAD::NearEqual;
    00780         bool ok  = true;
    00781         // first the objective function gradient
    00782         double checkObjPartial0 = 2*x0;
    00783         //ok &= NearEqual( *(objGrad + 0), checkObjPartial0, 1e-10, 1e-10);
    00784         ok = ( fabs(checkObjPartial0 - *(objGrad + 0) )/(fabs( checkObjPartial0) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00785         double checkObjPartial1 = 9;
    00786         //ok &= NearEqual( *(objGrad + 1), checkObjPartial1, 1e-10, 1e-10);
    00787         ok = ( fabs(checkObjPartial1 - *(objGrad + 1) )/(fabs( checkObjPartial1) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00788         double checkObjPartial2 = 0;
    00789         //ok &= NearEqual( *(objGrad + 2), checkObjPartial2, 1e-10, 1e-10);
    00790         ok = ( fabs(checkObjPartial2 - *(objGrad + 2) )/(fabs( checkObjPartial2) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00791         // get the constrating gradient
    00792         // row 0 gradient -- there are nonzero partials for variables 1 and 2
    00793         double checkCon0Partial1 = 1.37 + 5.0;
    00794         //ok &= NearEqual( *(sparseJac->values + 0), checkCon0Partial1, 1e-10, 1e-10);
    00795         ok = ( fabs(checkCon0Partial1 - *(sparseJac->values + 0) )/(fabs( checkCon0Partial1) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00796         double checkCon0Partial3 = 2.;
    00797         //ok &= NearEqual( *(sparseJac->values + 1), checkCon0Partial3, 1e-10, 1e-10);
    00798         ok = ( fabs(checkCon0Partial3 - *(sparseJac->values + 1) )/(fabs( checkCon0Partial3) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00799         // row 1 gradient -- there are nonzero partials for variables 0 and 2
    00800         double checkCon1Partial2 = 7;
    00801         //ok &= NearEqual( *(sparseJac->values + 2), checkCon1Partial2, 1e-10, 1e-10);
    00802         ok = ( fabs(checkCon1Partial2 - *(sparseJac->values + 2) )/(fabs( checkCon1Partial2) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00803         double checkCon1Partial0 = 1./x0;
    00804         //ok &= NearEqual( *(sparseJac->values + 3), checkCon1Partial0, 1e-10, 1e-10);
    00805         ok = ( fabs(checkCon1Partial0 - *(sparseJac->values + 3) )/(fabs( checkCon1Partial0) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00806         double checkCon1Partial3 = 1./x3;
    00807         //ok &= NearEqual( *(sparseJac->values + 4), checkCon1Partial3, 1e-10, 1e-10);
    00808         ok = ( fabs(checkCon1Partial3 - *(sparseJac->values + 4) )/(fabs( checkCon1Partial3) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
    00809         return ok;
    00810 }//CheckGradientValues
    00811 //
    00812 bool CheckHessianUpper( SparseHessianMatrix *sparseHessian ,
    00813         double x0, double x1, double x2, double x3, double z0, double z1, double w ){
    00814         using CppAD::NearEqual;
    00815         bool ok  = true;
    00816         int hessValuesIdx = 0;
    00817         //assert( sparseHessian->hessDimension = n * (n + 1) /2)
    00818         // L_00 = 2 * w - z1 / ( x0 * x0 )
    00819         double check = 2. * w - z1 / (x0 * x0);
    00820         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), check, 1e-10, 1e-10);
    00821         if(ok == false) std::cout << "FAILED ONE" << std::endl;
    00822         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), 0., 1e-10, 1e-10);
    00823         if(ok == false) std::cout << "FAILED TWO" << std::endl;
    00824         // L_22 = - z1 / (x3 * x3)
    00825         check = - z1 / (x3 * x3);
    00826         ok &= NearEqual(*(sparseHessian->hessValues + hessValuesIdx++), check, 1e-10, 1e-1);
    00827         if(ok == false) std::cout << "FAILED THREE" << std::endl;
    00828         return ok;
    00829 }//CheckHessianUpper