00001
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
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
00110 using CppAD::vector;
00111 std::cout.precision(12);
00112
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
00123
00124 FileUtil *fileUtil = NULL;
00125 std::string osilFileName;
00126 std::string osil;
00127
00128 const char dirsep = CoinFindDirSeparator();
00129
00130 std::string dataDir;
00131 dataDir = dirsep == '/' ? "../data/" : "..\\data\\";
00132
00133 osilFileName = dataDir + "CppADTestLag.osil";
00134 std::cout << "osilFileName = " << osilFileName << std::endl;
00135 fileUtil = new FileUtil();
00136 osil = fileUtil->getFileAsString( &osilFileName[0]);
00137
00138
00139 OSiLReader *osilreader = NULL;
00140 OSInstance *osinstance = NULL;
00141
00142 try{
00143
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
00151 double *objVals = NULL;
00152
00153
00167
00168
00169 std::cout << "Initialize Nonlinear Structures" << std::endl;
00170 osinstance->initForAlgDiff( );
00171
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
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
00206
00207
00208
00209
00210 size_t n = varIndexMap.size();
00211
00212 size_t m = 3;
00213
00214 std::vector<double> x0( n);
00219 x0[0] = 1;
00220 x0[1] = 5;
00221 x0[2] = 5;
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
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
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
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
00280
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
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];
00325 double* z = new double[2];
00326 double* w = new double[1];
00327 x[ 0] = 1;
00328 x[ 1] = 5;
00329 x[ 2] = 10;
00330 x[ 3] = 5;
00331 z[ 0] = 2;
00332 z[ 1] = 1;
00333 w[ 0] = 1;
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
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
00405
00406
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
00416
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
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
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
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
00489
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
00555 z[ 0] = 0;
00556 z[ 1] = 1;
00557 w[ 0] = 0;
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
00567 x[ 0] = 1;
00568
00569
00570
00571
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
00582 std::cout << "Start Taping" << std::endl;
00583 CppAD::Independent( X);
00584
00585
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
00590 CppAD::ADFun<double> f(X, Y);
00591 std::cout << "Stop Taping" << std::endl;
00592
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
00614
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
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
00630 gradVals = f.Forward(1, unit_col_vec);
00631 unit_col_vec[ index] = 0;
00632
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
00638 std::cout << "CALL f.Reverse -------" << std::endl;
00639 f.Reverse(2, lagMultipliers);
00640 std::cout << "FINISH CALL f.Reverse -------" << std::endl;
00641 }
00642
00643
00644 delete osilreader;
00645 osilreader = NULL;
00646 std::cout << "OSILREADER DELETED" << std::endl;
00647
00648
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
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
00667 CppAD::Independent(x);
00668
00669 size_t m = 1;
00670 CppADvector< AD<double> > y(m);
00671 y[0] = std::pow(x0, x1);
00672
00673 CppAD::ADFun<double> f(x, y);
00674
00675 double check = std::pow(x0, x1);
00676
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
00685 ok = ( fabs(check - dy[0])/(fabs( check) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
00686 }
00687
00688 {
00689
00690
00691
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
00698 CppAD::Independent(X);
00699
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
00705 CppAD::ADFun<double> f(X, Y);
00706
00707
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
00716 f.ForSparseJac(n, r);
00720
00721
00722
00723
00724
00725
00726
00727
00731
00732 std::vector<bool> e( m);
00733
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
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 }
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
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
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
00772 ok = ( fabs(checkCon1 - *(conVals + 1) )/(fabs( checkCon1) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
00773 return ok;
00774 }
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
00782 double checkObjPartial0 = 2*x0;
00783
00784 ok = ( fabs(checkObjPartial0 - *(objGrad + 0) )/(fabs( checkObjPartial0) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
00785 double checkObjPartial1 = 9;
00786
00787 ok = ( fabs(checkObjPartial1 - *(objGrad + 1) )/(fabs( checkObjPartial1) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
00788 double checkObjPartial2 = 0;
00789
00790 ok = ( fabs(checkObjPartial2 - *(objGrad + 2) )/(fabs( checkObjPartial2) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
00791
00792
00793 double checkCon0Partial1 = 1.37 + 5.0;
00794
00795 ok = ( fabs(checkCon0Partial1 - *(sparseJac->values + 0) )/(fabs( checkCon0Partial1) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
00796 double checkCon0Partial3 = 2.;
00797
00798 ok = ( fabs(checkCon0Partial3 - *(sparseJac->values + 1) )/(fabs( checkCon0Partial3) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
00799
00800 double checkCon1Partial2 = 7;
00801
00802 ok = ( fabs(checkCon1Partial2 - *(sparseJac->values + 2) )/(fabs( checkCon1Partial2) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
00803 double checkCon1Partial0 = 1./x0;
00804
00805 ok = ( fabs(checkCon1Partial0 - *(sparseJac->values + 3) )/(fabs( checkCon1Partial0) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
00806 double checkCon1Partial3 = 1./x3;
00807
00808 ok = ( fabs(checkCon1Partial3 - *(sparseJac->values + 4) )/(fabs( checkCon1Partial3) + OS_NEAR_EQUAL) <= OS_NEAR_EQUAL) ? true : false;
00809 return ok;
00810 }
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
00818
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
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 }