/home/coin/SVN-release/OS-1.1.1/OS/src/OSCommonInterfaces/OSInstance.cpp

Go to the documentation of this file.
00001 
00018 #include "OSInstance.h"
00019 #include "OSMathUtil.h"
00020 #include "OSCommonUtil.h"
00021 #include "OSErrorClass.h"
00022 #include "OSParameters.h"
00023 #include<iostream>  
00024 #include<sstream>
00025 
00026 //#define DEBUG
00027  
00028 using namespace std;
00029 using CppAD::NearEqual;
00030 
00031 OSInstance::OSInstance(): 
00032         m_sInstanceName(""),
00033         m_sInstanceSource(""),  
00034         m_sInstanceDescription(""),
00035         m_bProcessVariables(false),
00036         m_iVariableNumber(-1),
00037         m_iNumberOfIntegerVariables( 0),
00038         m_iNumberOfBinaryVariables( 0),
00039         m_iNumberOfQuadraticRowIndexes( 0),
00040         m_bQuadraticRowIndexesProcessed( false),
00041         m_miQuadRowIndexes( NULL),
00042         m_iNumberOfNonlinearExpressionTreeIndexes( 0),
00043         m_bNonlinearExpressionTreeIndexesProcessed( false),
00044         m_miNonlinearExpressionTreeIndexes( NULL),
00045         m_iNumberOfNonlinearExpressionTreeModIndexes( 0),
00046         m_bNonlinearExpressionTreeModIndexesProcessed( false),
00047         m_miNonlinearExpressionTreeModIndexes( NULL),           
00048         m_msVariableNames(NULL),
00049         m_mdVariableInitialValues(NULL),
00050         m_msVariableInitialStringValues(NULL),
00051         m_mcVariableTypes(NULL),
00052         m_mdVariableLowerBounds(NULL),
00053         m_mdVariableUpperBounds(NULL),
00054         m_bProcessObjectives(false),
00055         m_iObjectiveNumber(-1),
00056         m_iObjectiveNumberNonlinear( 0),
00057         m_msObjectiveNames(NULL),
00058         m_msMaxOrMins(NULL),
00059         m_miNumberOfObjCoef(NULL),
00060         m_mdObjectiveConstants(NULL),
00061         m_mdObjectiveWeights(NULL),
00062         m_mObjectiveCoefficients(NULL),
00063         m_bGetDenseObjectives(false),
00064         m_mmdDenseObjectiveCoefficients(NULL),
00065         m_bProcessConstraints(false),
00066         m_iConstraintNumber(-1),
00067         m_iConstraintNumberNonlinear( 0),       
00068         m_msConstraintNames(NULL),
00069         m_mdConstraintLowerBounds(NULL),
00070         m_mdConstraintUpperBounds(NULL),
00071         m_mdConstraintConstants( NULL),
00072         m_mcConstraintTypes(NULL),
00073         m_bProcessLinearConstraintCoefficients(false),  
00074         m_iLinearConstraintCoefficientNumber(-1),
00075         m_bColumnMajor(true),
00076         m_binitForAlgDiff( false),      
00077         m_linearConstraintCoefficientsInColumnMajor(NULL),
00078         m_linearConstraintCoefficientsInRowMajor(NULL), 
00079         m_bProcessQuadraticTerms(false),
00080         m_iQuadraticTermNumber(-1),
00081         m_mdConstraintFunctionValues( NULL),
00082         m_mdObjectiveFunctionValues( NULL),
00083         m_iJacValueSize( 0),
00084         m_miJacStart( NULL),
00085         m_miJacIndex( NULL),
00086         m_mdJacValue( NULL),
00087         m_miJacNumConTerms( NULL),
00088         m_sparseJacMatrix( NULL),       
00089         m_iHighestTaylorCoeffOrder(-1), 
00090         m_quadraticTerms( NULL),        
00091         m_bQTermsAdded( false), 
00092         m_iNumberOfNonlinearVariables( 0),
00093         m_bProcessNonlinearExpressions( false),
00094         m_iNonlinearExpressionNumber( -1),              
00095         m_miNonlinearExpressionIndexes( NULL),
00096         m_bProcessExpressionTrees( false),
00097         m_bProcessExpressionTreesMod( false),
00098         m_LagrangianExpTree(NULL),
00099         m_bLagrangianExpTreeCreated( false),
00100         m_LagrangianSparseHessian( NULL),
00101         m_bLagrangianSparseHessianCreated( false),
00102         m_miNonLinearVarsReverseMap( NULL),
00103         m_bAllNonlinearVariablesIndex( false),
00104         m_bCppADFunIsCreated( false),
00105         m_bCppADTapesBuilt( false),
00106         m_bCppADMustReTape( false),
00107         m_bDuplicateExpressionTreesMap( false),
00108         m_bNonLinearStructuresInitialized( false),
00109         m_bSparseJacobianCalculated( false),
00110         m_iHighestOrderEvaluated( -1),
00111         m_mmdObjGradient( NULL),
00112         m_bProcessTimeDomain( false),
00113         m_bProcessTimeStages( false),
00114         m_bProcessTimeInterval( false),
00115         m_bFiniteTimeStages( false),
00116         m_iNumberOfTimeStages(-1),
00117         bUseExpTreeForFunEval( false)
00118 
00119 {    
00120         #ifdef DEBUG
00121         cout << "Inside OSInstance Constructor" << endl;
00122         #endif
00123         this->instanceHeader = new InstanceHeader();
00124         this->instanceData = new InstanceData();
00125 }  
00126 
00127 OSInstance::~OSInstance(){
00128         #ifdef DEBUG
00129         cout << "OSInstance Destructor Called" << endl;
00130         #endif
00131         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
00132         // delete  the temporary arrays
00133         delete[] m_msVariableNames;
00134         m_msVariableNames = NULL;
00135         delete[] m_mdVariableInitialValues;
00136         m_mdVariableInitialValues = NULL ;
00137         delete[] m_msVariableInitialStringValues;
00138         m_msVariableInitialStringValues = NULL;
00139         delete[] m_mcVariableTypes;
00140         m_mcVariableTypes = NULL;
00153         delete[] m_miNonLinearVarsReverseMap;
00154         m_miNonLinearVarsReverseMap = NULL;
00155         int i;
00156         //if(instanceData->objectives->numberOfObjectives > 0 && m_mObjectiveCoefficients != NULL){
00157         if(m_bProcessObjectives == true ){
00158                 for(i = 0; i < instanceData->objectives->numberOfObjectives; i++){
00159                         #ifdef DEBUG
00160                         std::cout <<  "Delete m_mObjectiveCoefficients[i]" << std::endl;
00161                         #endif
00162                         delete m_mObjectiveCoefficients[i];
00163                         m_mObjectiveCoefficients[i] = NULL;
00164                 }
00165                 #ifdef DEBUG
00166                 std::cout <<  "Delete m_msObjectiveNames" << std::endl;
00167                 std::cout <<  "Delete m_msMaxOrMins" << std::endl;
00168                 std::cout <<  "Delete m_miNumberOfObjCoef" << std::endl;
00169                 std::cout <<  "Delete m_mdObjectiveConstants" << std::endl;
00170                 std::cout <<  "Delete m_mdObjectiveWeights" << std::endl;
00171                 #endif          
00172                 delete[] m_msObjectiveNames;
00173                 m_msObjectiveNames = NULL;
00174                 delete[] m_msMaxOrMins;
00175                 m_msMaxOrMins = NULL;
00176                 delete[] m_miNumberOfObjCoef;
00177                 m_miNumberOfObjCoef = NULL;
00178                 delete[] m_mdObjectiveConstants; 
00179                 m_mdObjectiveConstants = NULL;
00180                 delete[] m_mdObjectiveWeights;
00181                 m_mdObjectiveWeights = NULL;
00182                 delete[] m_mObjectiveCoefficients;
00183                 m_mObjectiveCoefficients = NULL;
00184         }
00185         if(instanceData->objectives->numberOfObjectives > 0 && m_mmdDenseObjectiveCoefficients != NULL){
00186                 for(i = 0; i < instanceData->objectives->numberOfObjectives; i++){
00187                         //delete m_mmdDenseObjectiveCoefficients[i];
00188                         #ifdef DEBUG
00189                         std::cout <<  "delete m_mmdDenseObjectiveCoefficients[i]" << std::endl;
00190                         #endif
00191                     delete[] m_mmdDenseObjectiveCoefficients[i];
00192                         m_mmdDenseObjectiveCoefficients[i] = NULL;
00193                 }
00194                 delete[] m_mmdDenseObjectiveCoefficients;
00195                 m_mmdDenseObjectiveCoefficients = NULL;
00196         }
00197 
00198         if( (m_binitForAlgDiff == true)  ){     
00199                 if(instanceData->objectives->numberOfObjectives > 0 && m_mmdObjGradient != NULL){
00200                         
00201                         #ifdef DEBUG
00202                         std::cout <<  "The number of objectives =  " << instanceData->objectives->numberOfObjectives << std::endl;
00203                         #endif
00204                         for(i = 0; i < instanceData->objectives->numberOfObjectives; i++){
00205                                 #ifdef DEBUG
00206                                 std::cout << "deleting Objective function gradient " << i << std::endl;
00207                                 #endif
00208                                 delete[] m_mmdObjGradient[i];
00209 
00210                                 m_mmdObjGradient[i] = NULL;
00211                         }
00212                         delete[] m_mmdObjGradient;
00213                         m_mmdObjGradient = NULL;
00214                 }
00215         }
00216 
00217         if(m_bProcessLinearConstraintCoefficients == true && m_bColumnMajor == true) delete m_linearConstraintCoefficientsInColumnMajor;
00218         if(m_bProcessLinearConstraintCoefficients == true && m_bColumnMajor == false) delete m_linearConstraintCoefficientsInRowMajor;
00219         //if(m_linearConstraintCoefficientsInRowMajor != NULL) delete m_linearConstraintCoefficientsInRowMajor;
00220         //if(m_linearConstraintCoefficientsInColumnMajor != NULL) delete m_linearConstraintCoefficientsInColumnMajor;
00221         delete[] m_msConstraintNames;
00222         m_msConstraintNames = NULL;
00223         delete[] m_mcConstraintTypes;
00224         m_mcConstraintTypes = NULL;
00225         delete[]  m_mdConstraintConstants;
00226          m_mdConstraintConstants = NULL;
00227         delete[] m_mdConstraintLowerBounds;
00228         m_mdConstraintLowerBounds = NULL;
00229         delete[] m_mdConstraintUpperBounds;
00230         m_mdConstraintUpperBounds = NULL;
00231         delete[] m_mdVariableLowerBounds;
00232         m_mdVariableLowerBounds = NULL;
00233         delete[] m_mdVariableUpperBounds;
00234         m_mdVariableUpperBounds = NULL;
00235         //std::cout << "Do garbage collection for the nonlinear API" << std::endl;
00236         // garbage collection for the gradient
00237         if(m_bNonLinearStructuresInitialized == true ){
00238                 delete[] m_mdObjectiveFunctionValues;
00239                 m_mdObjectiveFunctionValues = NULL;     
00240                 delete[] m_mdConstraintFunctionValues;
00241                 m_mdConstraintFunctionValues = NULL;
00242         }
00243         if(m_bSparseJacobianCalculated == true){
00244                 delete[] m_miJacStart;
00245                 m_miJacStart = NULL;
00246                 delete[] m_miJacIndex;
00247                 m_miJacIndex = NULL;
00248                 delete[] m_mdJacValue;
00249                 m_mdJacValue = NULL;
00250                 delete[] m_miJacNumConTerms;
00251                 m_miJacNumConTerms = NULL;
00252         }
00253         if( m_bLagrangianExpTreeCreated == true){
00254                 delete m_LagrangianExpTree;
00255                 m_LagrangianExpTree = NULL;
00256         }
00257         if( m_bLagrangianSparseHessianCreated == true){
00258                 delete m_LagrangianSparseHessian;
00259                 m_LagrangianSparseHessian = NULL;
00260         }
00261         if( m_bSparseJacobianCalculated == true){
00262                 delete m_sparseJacMatrix;
00263                 m_sparseJacMatrix = NULL;
00264         }
00265         if( (instanceData->quadraticCoefficients->qTerm != NULL) && (m_bProcessQuadraticTerms == true) ){
00266                 delete m_quadraticTerms;
00267                 m_quadraticTerms = NULL;
00268         }
00269         if( (instanceData->quadraticCoefficients->qTerm != NULL)  && (m_bQuadraticRowIndexesProcessed == true) ){
00270                 delete[] m_miQuadRowIndexes;
00271                 m_miQuadRowIndexes = NULL;
00272         }
00273         //
00274         // delete the new expression trees that got created
00275         //if( m_bLagrangianExpTreeCreated == false  ||  m_bLagrangianExpTreeCreated == true){
00276         if( (m_bProcessExpressionTrees == true) && (m_bDuplicateExpressionTreesMap == false)  ) {
00277                 for(posMapExpTree = m_mapExpressionTrees.begin(); posMapExpTree != m_mapExpressionTrees.end(); ++posMapExpTree){
00278                         std::cout << "Deleting an expression tree from the map for row  " << posMapExpTree->first  << std::endl;
00279                         delete m_mapExpressionTrees[ posMapExpTree->first ];
00280                 }
00281         }
00282         if( m_bDuplicateExpressionTreesMap == true)   {
00283                 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){          
00284                         #ifdef DEBUG
00285                                 std::cout << "Deleting an expression tree from m_mapExpressionTreesMod" << std::endl;
00286                         #endif
00287                         delete m_mapExpressionTreesMod[ posMapExpTree->first ];
00288                 }
00289         }
00290         //}
00292         if( (m_bNonlinearExpressionTreeIndexesProcessed == true) && (m_mapExpressionTrees.size() > 0) ){
00293                 std::cout << "Deleting  m_miNonlinearExpressionTreeIndexes" << std::endl;
00294                 delete[] m_miNonlinearExpressionTreeIndexes;
00295                 std::cout << "Done Deleting  m_miNonlinearExpressionTreeIndexes" << std::endl;
00296                 m_miNonlinearExpressionTreeIndexes = NULL;
00297         }
00298         if( (m_bNonlinearExpressionTreeModIndexesProcessed == true) && (m_mapExpressionTreesMod.size() > 0) ){
00299                 std::cout << "Deleting  m_miNonlinearExpressionTreeModIndexes" << std::endl;
00300                 delete[] m_miNonlinearExpressionTreeModIndexes;
00301                 std::cout << "Done Deleting  m_miNonlinearExpressionTreeModIndexes" << std::endl;
00302                 m_miNonlinearExpressionTreeModIndexes = NULL;
00303         }
00304         if(m_bCppADFunIsCreated == true){
00305                 delete Fad;
00306                 Fad = NULL;
00307         }
00308 //      if( (instanceData->timeDomain->stages->stage != NULL) && (m_bProcessTimeStages == true) ){
00309 //              delete m_Stages;
00310 //              m_Stages = NULL;
00311 //      }
00312         
00313         // delete the two children of OSInstance
00314         //delete instanceHeader object
00315         delete instanceHeader;
00316         instanceHeader = NULL;
00317         //delete instanceData object
00318         delete instanceData;
00319         instanceData = NULL;
00320 }//OSInstance Destructor
00321 
00322 InstanceHeader::InstanceHeader():
00323         description(""),
00324         name(""),
00325         source("")
00326 
00327 { 
00328         #ifdef DEBUG
00329         cout << "Inside the InstanceHeader Constructor" << endl;
00330         #endif
00331 } 
00332 
00333 
00334 InstanceHeader::~InstanceHeader(){
00335         #ifdef DEBUG  
00336         cout << "Inside the InstanceHeader Destructor" << endl;
00337         #endif
00338 } 
00339 
00340 Variable::Variable():
00341         lb(0.0),
00342         ub(OSDBL_MAX),
00343         init(OSNAN), 
00344         type('C'), 
00345         name(""),
00346         initString("")
00347 {  
00348         #ifdef DEBUG
00349         cout << "Inside the Variable Constructor" << endl;
00350         #endif 
00351 } 
00352 
00353 Variable::~Variable(){  
00354         #ifdef DEBUG
00355         cout << "Inside the Variable Destructor" << endl; 
00356         #endif
00357 } 
00358 
00359 Variables::Variables(){  
00360         #ifdef DEBUG
00361         cout << "Inside the Variables Constructor" << endl; 
00362         #endif 
00363         numberOfVariables = 0;
00364         var = NULL; 
00365 }
00366 
00367 Variables::~Variables(){ 
00368         #ifdef DEBUG 
00369         cout << "Inside the Variables Destructor" << endl;
00370         #endif
00371         int i;
00372         if(numberOfVariables > 0 && var != NULL){
00373                 for(i = 0; i < numberOfVariables; i++){
00374                         #ifdef DEBUG 
00375                         cout << "Deleting var[ i]" << endl;
00376                         #endif
00377                         delete var[i];
00378                         var[i] = NULL;
00379                 }
00380         }
00381         delete[] var;
00382         var = NULL; 
00383 }  
00384 
00385 ObjCoef::ObjCoef():
00386         idx(-1),
00387         value(0.0)  
00388 {  
00389         #ifdef DEBUG
00390         cout << "Inside the Coef Constructor" << endl;
00391         #endif 
00392 }
00393 
00394 ObjCoef::~ObjCoef(){ 
00395         #ifdef DEBUG
00396         cout << "Inside the ObjCoef Destructor" << endl;  
00397         #endif
00398 }
00399 
00400 Objective::Objective():
00401         name("") ,
00402         maxOrMin("min"),
00403         constant(0.0),
00404         weight(1.0),
00405         numberOfObjCoef(0),
00406         coef(NULL)
00407 { 
00408  
00409         #ifdef DEBUG
00410         cout << "Inside the Objective Constructor" << endl;
00411         #endif
00412 }
00413 
00414 Objective::~Objective(){
00415         #ifdef DEBUG  
00416         cout << "Inside the Objective Destructor" << endl;
00417         #endif
00418         int i;
00419         if(numberOfObjCoef > 0 && coef != NULL){
00420                 for(i = 0; i < numberOfObjCoef; i++){
00421                         delete coef[i];
00422                         coef[i] = NULL;
00423                 }
00424         }
00425         delete[] coef;
00426         coef = NULL;
00427 }  
00428 
00429 Objectives::Objectives()
00430 {  
00431         #ifdef DEBUG
00432         cout << "Inside the Objectives Constructor" << endl; 
00433         #endif
00434         numberOfObjectives = 0;
00435         obj = NULL;
00436 } 
00437 
00438 Objectives::~Objectives(){ 
00439         #ifdef DEBUG 
00440         cout << "Inside the Objectives Destructor" << endl;
00441         #endif
00442         int i;
00443         if(numberOfObjectives > 0 && obj != NULL){
00444                 for(i = 0; i < numberOfObjectives; i++){
00445                         delete obj[i];
00446                         obj[i] = NULL;
00447                 }
00448         }
00449         delete[] obj;
00450         obj = NULL;
00451 }
00452 
00453 Constraint::Constraint():
00454         name(""),
00455         constant(0.0),
00456         lb(-OSDBL_MAX),
00457         ub(OSDBL_MAX)
00458 
00459 {
00460         #ifdef DEBUG  
00461         cout << "Inside the Constraint Constructor" << endl;
00462         #endif
00463 } 
00464 
00465 Constraint::~Constraint(){  
00466         #ifdef DEBUG
00467         cout << "Inside the Constraint Destructor" << endl;
00468         #endif
00469 } 
00470 
00471 Constraints::Constraints():
00472         numberOfConstraints(0),
00473         con(NULL)
00474 {
00475         #ifdef DEBUG
00476         cout << "Inside the Constraints Constructor" << endl;
00477         #endif
00478 } 
00479 
00480 Constraints::~Constraints(){  
00481         #ifdef DEBUG
00482         cout << "Inside the Constraints Destructor" << endl;
00483         #endif
00484         int i;
00485         if(numberOfConstraints > 0 && con != NULL){
00486                 for( i = 0; i < numberOfConstraints; i++){
00487                         delete con[i];
00488                         con[i] = NULL;
00489                 }
00490         }
00491         delete[] con;
00492         con = NULL;
00493 } 
00494 
00495 
00496 
00497 LinearConstraintCoefficients::LinearConstraintCoefficients():
00498         numberOfValues(0) ,
00499         iNumberOfStartElements( 0)
00500 { 
00501         #ifdef DEBUG 
00502         cout << "Inside the LinearConstraintCoefficients Constructor" << endl; 
00503         #endif
00504         start = new IntVector();
00505         rowIdx = new IntVector();
00506         colIdx = new IntVector();
00507         value = new DoubleVector();
00508 
00509 } 
00510 
00511 
00512 LinearConstraintCoefficients::~LinearConstraintCoefficients(){  
00513         #ifdef DEBUG
00514         cout << "Inside the LinearConstraintCoefficients Destructor" << endl; 
00515         #endif
00516         delete start;
00517         start = NULL;
00518         delete rowIdx;
00519         rowIdx = NULL;
00520         delete colIdx;
00521         colIdx = NULL;
00522         delete value;
00523         value = NULL;
00524 }
00525 
00526 QuadraticTerm::QuadraticTerm():
00527 
00528         idx(0),   
00529         idxOne(-1),
00530         idxTwo(-1),
00531         coef(0.0)
00532 
00533 {
00534         #ifdef DEBUG  
00535         cout << "Inside the QuadraticTerm Constructor" << endl;
00536         #endif
00537 } 
00538 
00539 
00540 QuadraticTerm::~QuadraticTerm(){  
00541         #ifdef DEBUG
00542         cout << "Inside the QuadraticTerm Destructor" << endl;
00543         #endif
00544 }
00545 
00546 
00547 
00548 QuadraticCoefficients::QuadraticCoefficients():
00549         numberOfQuadraticTerms(0),
00550         qTerm(NULL)
00551 { 
00552         #ifdef DEBUG 
00553         cout << "Inside the QuadraticCoefficients Constructor" << endl;
00554         #endif
00555 }//end QuadraticCoefficients() 
00556 
00557 
00558 QuadraticCoefficients::~QuadraticCoefficients(){
00559         #ifdef DEBUG  
00560         cout << "Inside the QuadraticCoefficients Destructor" << endl;
00561         #endif
00562         int i;
00563         if(numberOfQuadraticTerms > 0 && qTerm != NULL){
00564                 for( i = 0; i < numberOfQuadraticTerms; i++){
00565                         delete qTerm[i];
00566                         qTerm[i] = NULL;
00567                 }
00568         }
00569         delete[] qTerm;
00570         qTerm = NULL;  
00571 }//end ~QuadraticCoefficients()  
00572 
00573 
00574 Nl::Nl(){
00575         idx = 0;
00576         osExpressionTree = NULL;
00577         m_bDeleteExpressionTree = true;
00578 }//end Nl
00579  
00580  
00581 Nl::~Nl(){
00582         #ifdef DEBUG  
00583         cout << "Inside the Nl Destructor" << endl;
00584         #endif
00585         // don't delete the expression tree if we created a map of the expression
00586         // trees, otherwise we would destroy twice
00587         if( m_bDeleteExpressionTree == true){
00588                 delete osExpressionTree;
00589                 osExpressionTree = NULL;
00590         }
00591         
00592 }//end ~Nl
00593 
00594 
00595 
00596 NonlinearExpressions::NonlinearExpressions():
00597         numberOfNonlinearExpressions(0) ,
00598         nl(NULL)
00599 { 
00600         #ifdef DEBUG 
00601         cout << "Inside the NonlinearExpressions Constructor" << endl;
00602         #endif
00603 }//end NonlinearExpressions() 
00604 
00605 NonlinearExpressions::~NonlinearExpressions(){
00606         #ifdef DEBUG  
00607         cout << "Inside the NonlinearExpressions Destructor" << endl;
00608         cout << "NUMBER OF NONLINEAR EXPRESSIONS = " << numberOfNonlinearExpressions << endl;
00609         #endif
00610         int i;
00611         if(numberOfNonlinearExpressions > 0 && nl != NULL){
00612                 for( i = 0; i < numberOfNonlinearExpressions; i++){
00613                         #ifdef DEBUG  
00614                                 cout << "DESTROYING EXPRESSION " << nl[ i]->idx << endl;
00615                         #endif
00616                         delete nl[i];
00617                         nl[i] = NULL;
00618                 }
00619         }
00620         delete[] nl;
00621         nl = NULL;  
00622 }//end ~NonlinearExpressions()  
00623 
00624 
00625 TimeDomainStage::TimeDomainStage():
00626         name(""),
00627         nvar(0),
00628         ncon(0),
00629         nobj(0)
00630 { 
00631         #ifdef DEBUG 
00632         cout << "Inside the Stage Constructor" << endl;
00633         #endif
00634         variables   = NULL;
00635         constraints = NULL;
00636         objectives  = NULL;
00637 }//end TimeDomainStage() 
00638 
00639 
00640 TimeDomainStage::~TimeDomainStage(){
00641         #ifdef DEBUG  
00642         cout << "Inside the Stage Destructor" << endl;
00643         #endif
00644         if (variables != NULL)
00645         {       delete [] variables;
00646                 variables = NULL;
00647         }
00648         if (constraints != NULL)
00649         {       delete []  constraints;
00650                 constraints = NULL;
00651         }
00652         if (objectives != NULL)
00653         {       delete [] objectives;
00654                 objectives = NULL;
00655         }
00656 }//end ~TimeDomainStage()  
00657 
00658 
00659 TimeDomainStages::TimeDomainStages():
00660         numberOfStages(0),
00661         stage(NULL)
00662 {
00663         #ifdef DEBUG  
00664         cout << "Inside the Stages Constructor" << endl;
00665         #endif
00666 } 
00667 
00668 
00669 TimeDomainStages::~TimeDomainStages(){  
00670         #ifdef DEBUG
00671         cout << "Inside the Stages Destructor" << endl;
00672         #endif
00673         int i;
00674         if(numberOfStages > 0 && stage != NULL){
00675                 for( i = 0; i < numberOfStages; i++){
00676                         delete stage[i];
00677                         stage[i] = NULL;
00678                 }
00679         }
00680         delete[] stage;
00681         stage = NULL;  
00682 }
00683 
00684 TimeDomainInterval::TimeDomainInterval():
00685         intervalHorizon(0.0),
00686         intervalStart(0.0)
00687 {
00688         #ifdef DEBUG  
00689         cout << "Inside the Interval Constructor" << endl;
00690         #endif
00691 } 
00692 
00693 
00694 TimeDomainInterval::~TimeDomainInterval(){  
00695         #ifdef DEBUG
00696         cout << "Inside the Interval Destructor" << endl;
00697         #endif
00698 }
00699 
00700 TimeDomain::TimeDomain()
00701 {
00702         #ifdef DEBUG
00703         cout << "Inside the TimeDomain Constructor" << endl;
00704         #endif
00705         stages = NULL;
00706         interval = NULL;
00707 }
00708 
00709 TimeDomain::~TimeDomain()
00710 {  
00711         #ifdef DEBUG
00712         cout << "Inside the TimeDomain Destructor" << endl;
00713         #endif
00714         if (stages != NULL)
00715         {       delete stages;
00716                 stages = NULL;
00717         };
00718         if (interval != NULL)
00719         {       delete interval;
00720                 interval = NULL;
00721         };
00722 } 
00723 
00724 
00725 InstanceData::InstanceData(){ 
00726         #ifdef DEBUG 
00727         cout << "Inside the InstanceData Constructor" << endl;
00728         #endif 
00729         variables = new Variables();
00730         objectives = new Objectives();
00731         constraints = new Constraints();
00732         linearConstraintCoefficients = new LinearConstraintCoefficients();
00733         quadraticCoefficients = new QuadraticCoefficients();
00734         nonlinearExpressions = new NonlinearExpressions();
00735         timeDomain = NULL;
00736 } 
00737 
00738 InstanceData::~InstanceData(){  
00739         #ifdef DEBUG
00740         cout << "Inside the InstanceData Destructor" << endl; 
00741         #endif
00742         delete variables;
00743         variables = NULL;
00744         delete objectives;
00745         objectives = NULL;
00746         delete constraints;
00747         constraints = NULL;
00748         delete linearConstraintCoefficients;
00749         linearConstraintCoefficients = NULL;
00750         delete quadraticCoefficients;
00751         quadraticCoefficients = NULL;
00752         delete nonlinearExpressions;
00753         nonlinearExpressions = NULL;
00754         delete timeDomain;
00755         timeDomain = NULL;
00756 } 
00757 
00758 string OSInstance::getInstanceName(){
00759         if(  m_sInstanceName.length() <= 0){
00760                 m_sInstanceName = instanceHeader->name;
00761         }
00762         return m_sInstanceName;
00763 }//getInstanceName
00764 
00765 
00766 string OSInstance::getInstanceSource(){
00767         if( m_sInstanceSource.length() <= 0){
00768                 m_sInstanceSource = instanceHeader->source;
00769         }
00770         return m_sInstanceSource;
00771 }//getInstanceSource
00772 
00773 string OSInstance::getInstanceDescription(){
00774         if(m_sInstanceDescription.length() <= 0){
00775                 m_sInstanceDescription = instanceHeader->description;
00776         }
00777         return m_sInstanceDescription;
00778 }//getInstanceDescription
00779 
00780 int OSInstance::getVariableNumber(){
00781         if(m_iVariableNumber == -1){
00782                 m_iVariableNumber = instanceData->variables->numberOfVariables;
00783         }
00784         return m_iVariableNumber;
00785 }//getVariableNumber
00786 
00787 int OSInstance::getNumberOfNonlinearExpressions(){
00788         if(m_iNonlinearExpressionNumber == -1){
00789                 m_iNonlinearExpressionNumber = instanceData->nonlinearExpressions->numberOfNonlinearExpressions;
00790         }
00791         return m_iNonlinearExpressionNumber;
00792 }//getNumberOfNonlinearExpressions
00793 
00794 
00795 
00796 bool OSInstance::processVariables() {
00797         if(m_bProcessVariables) return true;
00798         m_bProcessVariables = true;
00799         string vartype ="CBIS";
00800         int i = 0;
00801         int n = instanceData->variables->numberOfVariables;
00802         try{
00803                 if(instanceData->variables->var == NULL) throw ErrorClass("no variables defined");
00804                 if(instanceData->variables->var[0]->name.length() > 0 || instanceData->variables->var[n-1]->name.length() > 0){
00805                         m_msVariableNames = new string[n];
00806                         for(i = 0; i < n; i++) m_msVariableNames[i] = instanceData->variables->var[i]->name;
00807                 } 
00808                         m_mdVariableInitialValues = new double[n];
00809                         for(i = 0; i < n; i++){
00810                                 if(CommonUtil::ISOSNAN(instanceData->variables->var[ 0]->init) == true ){                               
00811                                         m_mdVariableInitialValues[i] =  1.7171;
00812                                 }
00813                                 else{
00814                                         m_mdVariableInitialValues[i] = instanceData->variables->var[i]->init;
00815                                 }
00816                         }
00817                 //}
00818                 if((instanceData->variables->var[0]->initString.length() > 0)){
00819                         m_msVariableInitialStringValues = new string[n];
00820                         for(i = 0; i < n; i++) m_msVariableInitialStringValues[i] = instanceData->variables->var[i]->initString;
00821                 }
00822                 m_mcVariableTypes = new char[n];
00823                 m_mdVariableLowerBounds = new double[n];
00824                 m_mdVariableUpperBounds = new double[n];
00825                 for(i = 0; i < n; i++){
00826                         if(vartype.find(instanceData->variables->var[i]->type) == string::npos) throw ErrorClass("wrong variable type");
00827                         m_mcVariableTypes[i] = instanceData->variables->var[i]->type;
00828                         if(m_mcVariableTypes[i] == 'B') m_iNumberOfBinaryVariables++;
00829                         if(m_mcVariableTypes[i] == 'I') m_iNumberOfIntegerVariables++;
00830                         m_mdVariableLowerBounds[i] = instanceData->variables->var[i]->lb;
00831                         m_mdVariableUpperBounds[i] = instanceData->variables->var[i]->ub;
00832                 }
00833                 return true;
00834         } //end try
00835         catch(const ErrorClass& eclass){
00836                 throw ErrorClass( eclass.errormsg);
00837         } 
00838 }//processVariables
00839         
00840 string* OSInstance::getVariableNames() {
00841         processVariables();
00842         return m_msVariableNames;
00843 }//getVariableNames     
00844 
00845 double* OSInstance::getVariableInitialValues() {
00846         processVariables();
00847         return m_mdVariableInitialValues;
00848 }//getVariableInitialValues
00849 
00850 string* OSInstance::getVariableInitialStringValues() {
00851         processVariables();
00852         return m_msVariableInitialStringValues;
00853 }//getVariableInitialStringValues
00854 
00855 char* OSInstance::getVariableTypes() {
00856         processVariables();
00857         return m_mcVariableTypes;
00858 }//getVariableTypes
00859 
00860 int OSInstance::getNumberOfIntegerVariables() {
00861         processVariables();
00862         return m_iNumberOfIntegerVariables;
00863 }//getNumberOfIntegerVariables
00864 
00865 int OSInstance::getNumberOfBinaryVariables() {
00866         processVariables();
00867         return m_iNumberOfBinaryVariables;
00868 }//getNumberOfBinaryVariables
00869 
00870 double* OSInstance::getVariableLowerBounds() {
00871         processVariables();
00872         return m_mdVariableLowerBounds;
00873 }//getVariableLowerBounds
00874 
00875 double* OSInstance::getVariableUpperBounds() {
00876         processVariables();
00877         return m_mdVariableUpperBounds;
00878 }//getVariableUpperBounds
00879 
00880 int OSInstance::getObjectiveNumber(){
00881         if(m_iObjectiveNumber == -1){
00882                 m_iObjectiveNumber = instanceData->objectives->numberOfObjectives;
00883         }
00884         return m_iObjectiveNumber;
00885 }//getObjectiveNumber
00886 
00887 
00888 bool OSInstance::processObjectives() {
00889         if(m_bProcessObjectives) return true;
00890         m_bProcessObjectives = true;
00891         int i = 0;
00892         int j = 0;
00893         if(instanceData == NULL || instanceData->objectives == NULL || instanceData->objectives->obj == NULL || instanceData->objectives->numberOfObjectives == 0) return true;
00894         int n = instanceData->objectives->numberOfObjectives;
00895         try{
00896                 if(instanceData->objectives->obj[0]->name.length() > 0 || instanceData->objectives->obj[n-1]->name.length() > 0){
00897                         m_msObjectiveNames = new string[n];
00898                         for(i = 0; i < n; i++) m_msObjectiveNames[i] = instanceData->objectives->obj[i]->name;
00899                 }
00900                 m_msMaxOrMins = new string[n];
00901                 m_miNumberOfObjCoef = new int[n];
00902                 m_mdObjectiveConstants = new double[n];
00903                 m_mdObjectiveWeights = new double[n];
00904                 m_mObjectiveCoefficients = new SparseVector*[n];
00905                 for(i = 0; i < n; i++){
00906                         m_mObjectiveCoefficients[i] = new SparseVector(instanceData->objectives->obj[ j]->numberOfObjCoef);
00907                         //m_mObjectiveCoefficients[i]->bDeleteArrays=false;
00908                 }
00909                 
00910                 //for(i = 0; i < n; i++){
00911                 //      m_mObjectiveCoefficients[i] = new SparseVector();
00912                 //      m_mObjectiveCoefficients[i]->number = instanceData->objectives->obj[ j]->numberOfObjCoef;
00913                 //}
00914                 for(i = 0; i < n; i++){
00915                         if((instanceData->objectives->obj[i]->maxOrMin.compare("max") != 0) && (instanceData->objectives->obj[i]->maxOrMin.compare("min") != 0 )) throw ErrorClass("wrong objective maxOrMin");
00916                         m_msMaxOrMins[i] = instanceData->objectives->obj[i]->maxOrMin;
00917                         m_miNumberOfObjCoef[i] = instanceData->objectives->obj[i]->numberOfObjCoef;
00918                         m_mdObjectiveConstants[i] = instanceData->objectives->obj[i]->constant;
00919                         m_mdObjectiveWeights[i] = instanceData->objectives->obj[i]->weight;
00920                         if(instanceData->objectives->obj[i]->coef == NULL && m_miNumberOfObjCoef[i] != 0){
00921                                 throw ErrorClass("objective coefficient number inconsistent with objective coefficient array");
00922                         }
00923                         if(instanceData->objectives->obj[i]->coef != NULL){
00924                                 for(j = 0; j < m_mObjectiveCoefficients[i]->number; j++){
00925                                         m_mObjectiveCoefficients[i]->indexes[j] = instanceData->objectives->obj[i]->coef[j]->idx;
00926                                         m_mObjectiveCoefficients[i]->values[j] = instanceData->objectives->obj[i]->coef[j]->value;                      
00927                                 }
00928                         }
00929                 }               
00930                 return true;
00931         }
00932         catch(const ErrorClass& eclass){
00933                 throw ErrorClass( eclass.errormsg);
00934         }
00935 }//processObjectives
00936 
00937 string* OSInstance::getObjectiveNames() {
00938         processObjectives();
00939         return m_msObjectiveNames;
00940 }//getObjectiveNames
00941 
00942 string* OSInstance::getObjectiveMaxOrMins() {
00943         processObjectives();
00944         return m_msMaxOrMins;
00945 }//getObjectiveMaxOrMins
00946 
00947 int* OSInstance::getObjectiveCoefficientNumbers(){
00948                 processObjectives();
00949                 return m_miNumberOfObjCoef;
00950         }//getObjectiveCoefficientNumbers
00951 
00952 double* OSInstance::getObjectiveConstants() {
00953         processObjectives();
00954         return m_mdObjectiveConstants;
00955 }//getObjectiveConstants
00956 
00957 double* OSInstance::getObjectiveWeights() {
00958         processObjectives();
00959         return m_mdObjectiveWeights;
00960 }//getObjectiveWeights
00961 
00962 SparseVector** OSInstance::getObjectiveCoefficients() {
00963         processObjectives();
00964         return m_mObjectiveCoefficients;
00965 }//getObjectiveCoefficients
00966 
00967 
00968 double** OSInstance::getDenseObjectiveCoefficients() {
00969         if(m_bGetDenseObjectives) return m_mmdDenseObjectiveCoefficients;
00970         m_bGetDenseObjectives = true;
00971         if(instanceData->objectives->obj == NULL || instanceData->objectives->numberOfObjectives == 0) return m_mmdDenseObjectiveCoefficients;
00972         int m = instanceData->objectives->numberOfObjectives;
00973         int n = instanceData->variables->numberOfVariables;
00974         m_mmdDenseObjectiveCoefficients = new double*[m];
00975         int i, j, numobjcoef;
00976         SparseVector *sparsevec;
00977         for(i = 0; i < m; i++){
00978                 sparsevec = this->getObjectiveCoefficients()[i];
00979                 m_mmdDenseObjectiveCoefficients[ i] = new double[n];
00980                 for(j = 0; j < n; j++){
00981                         m_mmdDenseObjectiveCoefficients[ i][j] = 0.0;
00982                 }
00983                 sparsevec =  this->getObjectiveCoefficients()[i];
00984                 numobjcoef = sparsevec->number;
00985                 for(j = 0; j < numobjcoef; j++){
00986                         m_mmdDenseObjectiveCoefficients[i][ sparsevec->indexes[ j]]
00987                         = sparsevec->values[ j];
00988                 }
00989         }
00990         return m_mmdDenseObjectiveCoefficients;
00991 }//getDenseObjectiveCoefficients
00992 
00993 int OSInstance::getConstraintNumber(){
00994         if(m_iConstraintNumber == -1){
00995                 m_iConstraintNumber = instanceData->constraints->numberOfConstraints;
00996         }
00997         return m_iConstraintNumber;
00998 }//getConstraintNumber
00999 
01000 bool OSInstance::processConstraints() {
01001         if(m_bProcessConstraints) return true;
01002         m_bProcessConstraints = true;
01003         int i = 0;
01004         ostringstream outStr;
01005         if(instanceData == NULL || instanceData->constraints == NULL || instanceData->constraints->con == NULL || instanceData->constraints->numberOfConstraints == 0) return true;
01006         int n = instanceData->constraints->numberOfConstraints;
01007         try{
01008                 if(instanceData->constraints->con[0]->name.length() > 0 || instanceData->constraints->con[n-1]->name.length() > 0){
01009                         m_msConstraintNames = new string[n];
01010                         for(i = 0; i < n; i++) m_msConstraintNames[i] = instanceData->constraints->con[i]->name;
01011                 }
01012                 m_mdConstraintLowerBounds = new double[n];
01013                 m_mdConstraintUpperBounds = new double[n];
01014                 m_mdConstraintConstants = new double[n];
01015                 m_mcConstraintTypes = new char[n];
01016                 for(i = 0; i < n; i++){
01017                         m_mdConstraintLowerBounds[i] = instanceData->constraints->con[i]->lb;
01018                         m_mdConstraintUpperBounds[i] = instanceData->constraints->con[i]->ub;
01019                         m_mdConstraintConstants[i] = instanceData->constraints->con[i]->constant;
01020                         if(m_mdConstraintLowerBounds[i] == OSDBL_MAX || m_mdConstraintUpperBounds[i] == -OSDBL_MAX) {
01021                                 throw ErrorClass( outStr.str() );
01022                         }
01023                         else if(m_mdConstraintLowerBounds[i] > m_mdConstraintUpperBounds[i]) {
01024                                 outStr << "Constraint  " ;
01025                                 outStr << i;
01026                                 outStr << " is infeasible";
01027                                 throw ErrorClass( outStr.str());
01028                         }
01029                         else if(m_mdConstraintLowerBounds[i] == -OSDBL_MAX && m_mdConstraintUpperBounds[i] == OSDBL_MAX)
01030                                 m_mcConstraintTypes[i] = 'U';
01031                         else if(m_mdConstraintLowerBounds[i] == m_mdConstraintUpperBounds[i]) 
01032                                 m_mcConstraintTypes[i] = 'E';
01033                         else if(m_mdConstraintLowerBounds[i] == -OSDBL_MAX)
01034                                 m_mcConstraintTypes[i] = 'L';
01035                         else if(m_mdConstraintUpperBounds[i] == OSDBL_MAX)
01036                                 m_mcConstraintTypes[i] = 'G';
01037                         else m_mcConstraintTypes[i] = 'R';
01038                 }
01039                 return true;
01040         }
01041         catch(const ErrorClass& eclass){
01042                 throw ErrorClass( eclass.errormsg);
01043         }
01044 }//processConstraints
01045 
01046 
01047 string* OSInstance::getConstraintNames() {
01048         processConstraints();
01049         return m_msConstraintNames;
01050 }//getConstraintNames
01051 
01052 
01053 double* OSInstance::getConstraintLowerBounds() {
01054         processConstraints();
01055         return m_mdConstraintLowerBounds;
01056 }//getConstraintLowerBounds
01057 
01058 char* OSInstance::getConstraintTypes() {
01059         processConstraints();
01060         return m_mcConstraintTypes;
01061 }//getConstraintTypes
01062 
01063 double* OSInstance::getConstraintUpperBounds() {
01064         processConstraints();
01065         return m_mdConstraintUpperBounds;
01066 }//getConstraintUpperBounds
01067 
01068 int OSInstance::getLinearConstraintCoefficientNumber(){
01069         if(m_iLinearConstraintCoefficientNumber == -1){
01070                 m_iLinearConstraintCoefficientNumber = instanceData->linearConstraintCoefficients->numberOfValues;
01071         }
01072         return m_iLinearConstraintCoefficientNumber; 
01073 }//getLinearConstraintCoefficientNumber
01074 
01075 bool OSInstance::processLinearConstraintCoefficients() {
01076         if(m_bProcessLinearConstraintCoefficients) return true;
01077         m_bProcessLinearConstraintCoefficients = true;
01078         try{
01079                 int n = instanceData->linearConstraintCoefficients->numberOfValues;
01080                 //value array
01081                 if((instanceData->linearConstraintCoefficients->value == NULL ) || (n == 0) ) return true;
01082                 //index array
01083                 if((instanceData->linearConstraintCoefficients->colIdx != NULL && instanceData->linearConstraintCoefficients->colIdx->el != NULL) 
01084                 && (instanceData->linearConstraintCoefficients->rowIdx != NULL && instanceData->linearConstraintCoefficients->rowIdx->el != NULL))
01085                         throw ErrorClass("ambiguous linear constraint coefficient major");
01086                 else if(instanceData->linearConstraintCoefficients->value->el == NULL) return true;
01087                 else{
01088                         if(instanceData->linearConstraintCoefficients->rowIdx->el != NULL){
01089                                 m_bColumnMajor = true;
01090                                 m_linearConstraintCoefficientsInColumnMajor = new SparseMatrix();
01091                                 m_linearConstraintCoefficientsInColumnMajor->bDeleteArrays = false;
01092                                 m_linearConstraintCoefficientsInColumnMajor->isColumnMajor = true;
01093                                 m_linearConstraintCoefficientsInColumnMajor->valueSize = n;
01094                                 m_linearConstraintCoefficientsInColumnMajor->startSize = instanceData->variables->numberOfVariables + 1;
01095                         }
01096                         else{ 
01097                                 m_bColumnMajor = false; 
01098                                 m_linearConstraintCoefficientsInRowMajor = new SparseMatrix();
01099                                 m_linearConstraintCoefficientsInRowMajor->bDeleteArrays = false;
01100                                 m_linearConstraintCoefficientsInRowMajor->isColumnMajor = false;
01101                                 m_linearConstraintCoefficientsInRowMajor->valueSize = n;
01102                                 m_linearConstraintCoefficientsInRowMajor->startSize = instanceData->constraints->numberOfConstraints + 1;
01103                         }
01104                 }                       
01105                 if(m_bColumnMajor == true){
01106                         m_linearConstraintCoefficientsInColumnMajor->values = instanceData->linearConstraintCoefficients->value->el;
01107                         m_linearConstraintCoefficientsInColumnMajor->indexes = instanceData->linearConstraintCoefficients->rowIdx->el;
01108                         m_linearConstraintCoefficientsInColumnMajor->starts = instanceData->linearConstraintCoefficients->start->el;                    
01109                 }
01110                 else{
01111                         m_linearConstraintCoefficientsInRowMajor->values = instanceData->linearConstraintCoefficients->value->el;
01112                         m_linearConstraintCoefficientsInRowMajor->indexes = instanceData->linearConstraintCoefficients->colIdx->el;
01113                         m_linearConstraintCoefficientsInRowMajor->starts = instanceData->linearConstraintCoefficients->start->el;                                               
01114                 }
01115                 return true;
01116         }
01117         catch(const ErrorClass& eclass){
01118                 throw ErrorClass( eclass.errormsg);
01119         }
01120 }//processLinearConstraintCoefficients
01121 
01122 bool OSInstance::getLinearConstraintCoefficientMajor() {
01123         processLinearConstraintCoefficients();  
01124         return m_bColumnMajor;          
01125 }//getLinearConstraintCoefficientMajor
01126 
01127 SparseMatrix* OSInstance::getLinearConstraintCoefficientsInColumnMajor() {
01128         processLinearConstraintCoefficients();
01129         if(m_linearConstraintCoefficientsInColumnMajor != NULL) return m_linearConstraintCoefficientsInColumnMajor;
01130         if(!m_bColumnMajor){
01131                 if(m_linearConstraintCoefficientsInRowMajor == NULL) return NULL;
01132                 m_linearConstraintCoefficientsInColumnMajor = 
01133                         MathUtil::convertLinearConstraintCoefficientMatrixToTheOtherMajor(false,
01134                                           m_linearConstraintCoefficientsInRowMajor->startSize,
01135                                           m_linearConstraintCoefficientsInRowMajor->valueSize,
01136                                           m_linearConstraintCoefficientsInRowMajor->starts,
01137                                           m_linearConstraintCoefficientsInRowMajor->indexes,
01138                                           m_linearConstraintCoefficientsInRowMajor->values,
01139                                           getVariableNumber());
01140         }
01141         return m_linearConstraintCoefficientsInColumnMajor;             
01142 }//getLinearConstraintCoefficientsInColumnMajor
01143 
01144 SparseMatrix* OSInstance::getLinearConstraintCoefficientsInRowMajor() {
01145         processLinearConstraintCoefficients();
01146         if(m_linearConstraintCoefficientsInRowMajor != NULL) return m_linearConstraintCoefficientsInRowMajor;
01147         if(m_bColumnMajor){
01148                 if(m_linearConstraintCoefficientsInColumnMajor == NULL) return NULL;
01149                 m_linearConstraintCoefficientsInRowMajor = 
01150                 MathUtil::convertLinearConstraintCoefficientMatrixToTheOtherMajor(true,
01151                                   m_linearConstraintCoefficientsInColumnMajor->startSize,
01152                                   m_linearConstraintCoefficientsInColumnMajor->valueSize,
01153                                   m_linearConstraintCoefficientsInColumnMajor->starts,
01154                                   m_linearConstraintCoefficientsInColumnMajor->indexes,
01155                                   m_linearConstraintCoefficientsInColumnMajor->values,
01156                                   getConstraintNumber());
01157         }
01158         return m_linearConstraintCoefficientsInRowMajor; 
01159 }//getLinearConstraintCoefficientsInRowMajor
01160 
01161 
01162 int OSInstance::getNumberOfQuadraticTerms(){
01163         if(m_iQuadraticTermNumber == -1){
01164         // if m_iQuadraticTermNumber == -1 then the parser did not find any q terms so 
01165         // must new the object
01166                 if(instanceData->quadraticCoefficients == NULL)instanceData->quadraticCoefficients = new QuadraticCoefficients();
01167                 m_iQuadraticTermNumber = instanceData->quadraticCoefficients->numberOfQuadraticTerms;
01168         }
01169         return m_iQuadraticTermNumber;
01170 }//getNumberOfQuadraticTerms
01171 
01172 QuadraticTerms* OSInstance::getQuadraticTerms() {
01173         if(m_bProcessQuadraticTerms) return m_quadraticTerms;
01174         m_bProcessQuadraticTerms = true;
01175         if(instanceData->quadraticCoefficients->qTerm == 0) return 0;
01176         try{
01177                 int i = 0;
01178                 QuadraticCoefficients* quadraticCoefs = instanceData->quadraticCoefficients;
01179                 int n = quadraticCoefs->numberOfQuadraticTerms;
01180                 if(!quadraticCoefs->qTerm  && n != 0) 
01181                         throw ErrorClass("quadratic term number inconsistent with quadratic term array");               
01182                 m_quadraticTerms = new QuadraticTerms();
01183                 m_quadraticTerms->rowIndexes = new int[n];
01184                 m_quadraticTerms->varOneIndexes = new int[n];
01185                 m_quadraticTerms->varTwoIndexes = new int[n];
01186                 m_quadraticTerms->coefficients = new double[n];
01187                 for(i = 0; i < n; i++){
01188                         m_quadraticTerms->rowIndexes[i] = quadraticCoefs->qTerm[i]->idx;
01189                         m_quadraticTerms->varOneIndexes[i] = quadraticCoefs->qTerm[i]->idxOne;
01190                         m_quadraticTerms->varTwoIndexes[i] = quadraticCoefs->qTerm[i]->idxTwo;
01191                         m_quadraticTerms->coefficients[i] = quadraticCoefs->qTerm[i]->coef;
01192                 } 
01193                 return m_quadraticTerms;
01194         }
01195         catch(const ErrorClass& eclass){
01196                 throw ErrorClass( eclass.errormsg);
01197         } 
01198 }//getQuadraticTerms
01199 
01200 
01201 int* OSInstance::getQuadraticRowIndexes() {
01202         if(m_bQuadraticRowIndexesProcessed == true) return m_miQuadRowIndexes;
01203         m_bQuadraticRowIndexesProcessed = true;
01204         int n = getNumberOfQuadraticTerms();    
01205         if(n <= 0) return NULL;
01206         QuadraticTerms *qTerms = NULL;
01207         qTerms = getQuadraticTerms();
01208         std::map<int, int> foundIdx;
01209         std::map<int, int>::iterator pos;
01210         int i;
01211         try{
01212                 for(i = 0; i < n; i++){
01213                         // add the terms
01214                         foundIdx[ qTerms->rowIndexes[ i] ];      
01215                 }
01216                 // now put the term into an array
01217                 m_iNumberOfQuadraticRowIndexes = foundIdx.size();
01218                 m_miQuadRowIndexes = new int[ m_iNumberOfQuadraticRowIndexes ]  ;
01219                 i = 0;
01220                 for(pos = foundIdx.begin(); pos != foundIdx.end(); ++pos){
01221                         m_miQuadRowIndexes[ i++] = pos->first;  
01222                 }
01223                 foundIdx.clear();       
01224                 return m_miQuadRowIndexes;
01225         }
01226         catch(const ErrorClass& eclass){
01227                 throw ErrorClass( eclass.errormsg);
01228         } 
01229 }//getQuadraticRowIndexes
01230 
01231 
01232 int OSInstance::getNumberOfQuadraticRowIndexes() {
01233         if(m_bQuadraticRowIndexesProcessed == false) getQuadraticRowIndexes();
01234         return m_iNumberOfQuadraticRowIndexes;
01235 }//getNumberOfQuadraticRowIndexes
01236 
01237 int* OSInstance::getNonlinearExpressionTreeIndexes(){
01238         if(m_bNonlinearExpressionTreeIndexesProcessed == true) return m_miNonlinearExpressionTreeIndexes;
01239         m_bNonlinearExpressionTreeIndexesProcessed = true;
01240         std::map<int, OSExpressionTree*> expTrees;
01241         expTrees = getAllNonlinearExpressionTrees();    
01242         std::map<int, OSExpressionTree*>::iterator pos;
01243         try{
01244                 // now put the term into an array
01245                 m_iNumberOfNonlinearExpressionTreeIndexes = expTrees.size();
01246                 m_miNonlinearExpressionTreeIndexes = new int[ m_iNumberOfNonlinearExpressionTreeIndexes ]       ;
01247                 int i = 0;
01248                 for(pos = expTrees.begin(); pos != expTrees.end(); ++pos){
01249                         m_miNonlinearExpressionTreeIndexes[ i++] = pos->first;  
01250                 }
01251                 expTrees.clear();       
01252                 return m_miNonlinearExpressionTreeIndexes;
01253         }
01254         catch(const ErrorClass& eclass){
01255                 throw ErrorClass( eclass.errormsg);
01256         } 
01257 }//getNonlinearExpressionTreeIndexes
01258 
01259 int OSInstance::getNumberOfNonlinearExpressionTreeIndexes() {
01260         if(m_bNonlinearExpressionTreeIndexesProcessed == false) getNonlinearExpressionTreeIndexes();
01261         return m_iNumberOfNonlinearExpressionTreeIndexes;
01262 }//getNumberOfNonlinearExpressionTreeIndexes
01263 
01264 
01265 
01266 int* OSInstance::getNonlinearExpressionTreeModIndexes(){
01267         if(m_bNonlinearExpressionTreeModIndexesProcessed == true) return m_miNonlinearExpressionTreeModIndexes;
01268         m_bNonlinearExpressionTreeModIndexesProcessed = true;
01269         std::map<int, OSExpressionTree*> expTrees;
01270         expTrees = getAllNonlinearExpressionTreesMod(); 
01271         std::map<int, OSExpressionTree*>::iterator pos;
01272         try{
01273                 // now put the term into an array
01274                 m_iNumberOfNonlinearExpressionTreeModIndexes = expTrees.size();
01275                 m_miNonlinearExpressionTreeModIndexes = new int[ m_iNumberOfNonlinearExpressionTreeModIndexes ] ;
01276                 int i = 0;
01277                 for(pos = expTrees.begin(); pos != expTrees.end(); ++pos){
01278                         m_miNonlinearExpressionTreeModIndexes[ i++] = pos->first;       
01279                 }
01280                 expTrees.clear();       
01281                 return m_miNonlinearExpressionTreeModIndexes;
01282         }
01283         catch(const ErrorClass& eclass){
01284                 throw ErrorClass( eclass.errormsg);
01285         } 
01286 }//getNonlinearExpressionTreeModIndexes
01287 
01288 int OSInstance::getNumberOfNonlinearExpressionTreeModIndexes() {
01289         if(m_bNonlinearExpressionTreeModIndexesProcessed == false) getNonlinearExpressionTreeModIndexes();
01290         return m_iNumberOfNonlinearExpressionTreeModIndexes;
01291 }//getNumberOfNonlinearExpressionTreeModIndexes
01292 
01293 
01294 int OSInstance::getNumberOfNonlinearConstraints(){
01295         if( m_bProcessExpressionTrees == false ){
01296                 getAllNonlinearExpressionTrees();
01297                 return m_iConstraintNumberNonlinear;
01298         }
01299         else return m_iConstraintNumberNonlinear;
01300 }//getNumberOfNonlinearConstraints
01301 
01302 int OSInstance::getNumberOfNonlinearObjectives(){
01303         if( m_bProcessExpressionTrees == false ){
01304                 getAllNonlinearExpressionTrees();
01305                 return m_iObjectiveNumberNonlinear;
01306         }
01307         else return m_iObjectiveNumberNonlinear;
01308 }//getNumberOfNonlinearObjectivess
01309 
01310 
01311 OSExpressionTree* OSInstance::getNonlinearExpressionTree(int rowIdx){
01312         if( m_bProcessExpressionTrees == false ){
01313                 getAllNonlinearExpressionTrees();
01314                 return m_mapExpressionTrees[ rowIdx];
01315         } 
01316         else{
01317                 //if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()) return NULL;
01318                 //else return m_mapExpressionTrees[ rowIdx];
01319                 if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()) return m_mapExpressionTrees[ rowIdx];
01320                 else return NULL ;
01321                 // check to make sure rowIdx has a nonlinear term and is in the map
01323                 //std::map<int, OSExpressionTree*>::iterator pos;
01325                 //      if(pos->first == rowIdx)return m_mapExpressionTrees[ rowIdx];
01326                 //}
01327                 // if we are rowIdx has no nonlinar terms so return a null
01328                 //return NULL;
01329         }     
01330 }// getNonlinearExpressionTree for a specific index   
01331 
01332 
01333 OSExpressionTree* OSInstance::getNonlinearExpressionTreeMod(int rowIdx){
01334         if( m_bProcessExpressionTreesMod == false ){
01335                 getAllNonlinearExpressionTreesMod();
01336                 return m_mapExpressionTreesMod[ rowIdx];
01337         } 
01338         else{
01339                 //if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()) return NULL;
01340                 //else return m_mapExpressionTrees[ rowIdx];
01341                 if( m_mapExpressionTreesMod.find( rowIdx) != m_mapExpressionTreesMod.end()) return m_mapExpressionTreesMod[ rowIdx];
01342                 else return NULL ;
01343                 // check to make sure rowIdx has a nonlinear term and is in the map
01345                 //std::map<int, OSExpressionTree*>::iterator pos;
01347                 //      if(pos->first == rowIdx)return m_mapExpressionTrees[ rowIdx];
01348                 //}
01349                 // if we are rowIdx has no nonlinar terms so return a null
01350                 //return NULL;
01351         }     
01352 }// getNonlinearExpressionTreeMod for a specific index 
01353 
01354 
01355 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeInPostfix( int rowIdx){
01356         if( m_bProcessExpressionTrees == false ) getAllNonlinearExpressionTrees();
01357         std::vector<OSnLNode*> postfixVec;
01358         try{
01359                 if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()){
01360                         OSExpressionTree* expTree = getNonlinearExpressionTree( rowIdx);
01361                         postfixVec = expTree->m_treeRoot->getPostfixFromExpressionTree();
01362                         
01363                 }  
01364                 else{
01365                         throw ErrorClass("Error in getNonlinearExpressionTreeInPostfix, rowIdx not valid");
01366                 }
01367                 return postfixVec;      
01368         }
01369         catch(const ErrorClass& eclass){
01370                 throw ErrorClass( eclass.errormsg);
01371         } 
01372 }//getNonlinearExpressionTreeInPostfix
01373 
01374 
01375 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeModInPostfix( int rowIdx){
01376         if( m_bProcessExpressionTreesMod == false ) getAllNonlinearExpressionTreesMod();
01377         std::vector<OSnLNode*> postfixVec;
01378         try{
01379                 if( m_mapExpressionTreesMod.find( rowIdx) != m_mapExpressionTreesMod.end()){
01380                         OSExpressionTree* expTree = getNonlinearExpressionTreeMod( rowIdx);
01381                         postfixVec = expTree->m_treeRoot->getPostfixFromExpressionTree();
01382                         
01383                 }  
01384                 else{
01385                         throw ErrorClass("Error in getNonlinearExpressionTreeModInPostfix, rowIdx not valid");
01386                 }
01387                 return postfixVec;      
01388         }
01389         catch(const ErrorClass& eclass){
01390                 throw ErrorClass( eclass.errormsg);
01391         } 
01392 }//getNonlinearExpressionTreeModInPostfix
01393 
01394 
01395 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeInPrefix( int rowIdx){
01396         if( m_bProcessExpressionTrees == false ) getAllNonlinearExpressionTrees();
01397         std::vector<OSnLNode*> prefixVec;
01398         try{
01399                 if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()){
01400                         OSExpressionTree* expTree = getNonlinearExpressionTree( rowIdx);
01401                         prefixVec = expTree->m_treeRoot->getPrefixFromExpressionTree();
01402                         
01403                 }  
01404                 else{
01405                         throw ErrorClass("Error in getNonlinearExpressionTreeInPrefix, rowIdx not valid");
01406                 }
01407                 return prefixVec;       
01408         }
01409         catch(const ErrorClass& eclass){
01410                 throw ErrorClass( eclass.errormsg);
01411         } 
01412 }//getNonlinearExpressionTreeInPrefix
01413 
01414 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeModInPrefix( int rowIdx){
01415         if( m_bProcessExpressionTreesMod == false ) getAllNonlinearExpressionTreesMod();
01416         std::vector<OSnLNode*> prefixVec;
01417         try{
01418                 if( m_mapExpressionTreesMod.find( rowIdx) != m_mapExpressionTreesMod.end()){
01419                         OSExpressionTree* expTree = getNonlinearExpressionTreeMod( rowIdx);
01420                         prefixVec = expTree->m_treeRoot->getPrefixFromExpressionTree();
01421                         
01422                 }  
01423                 else{
01424                         throw ErrorClass("Error in getNonlinearExpressionTreeInPrefix, rowIdx not valid");
01425                 }
01426                 return prefixVec;       
01427         }
01428         catch(const ErrorClass& eclass){
01429                 throw ErrorClass( eclass.errormsg);
01430         } 
01431 }//getNonlinearExpressionTreeInPrefix
01432 
01433 std::map<int, OSExpressionTree*> OSInstance::getAllNonlinearExpressionTrees(){
01434         if(m_bProcessExpressionTrees == true) return m_mapExpressionTrees;
01435         std::map<int, int> foundIdx;
01436         std::map<int, int>::iterator pos;
01437         OSnLNodePlus *nlNodePlus;
01438         OSExpressionTree *expTree;  
01439         m_iObjectiveNumberNonlinear = 0;   
01440         m_iConstraintNumberNonlinear = 0;    
01441         int i;   
01442         // important -- tell the nl nodes not to destroy the OSExpression Objects
01443         if( instanceData->nonlinearExpressions->numberOfNonlinearExpressions > 0 && instanceData->nonlinearExpressions->nl != NULL){
01444                 for( i = 0; i < instanceData->nonlinearExpressions->numberOfNonlinearExpressions; i++){
01445                         instanceData->nonlinearExpressions->nl[i]->m_bDeleteExpressionTree = false;
01446                 }
01447         }
01448         int index;  
01449         // kipp -- what should we return if instanceData->nonlinearExpressions->numberOfNonlinearExpressions is zero
01450         for(i = 0; i < instanceData->nonlinearExpressions->numberOfNonlinearExpressions; i++){
01451                 index = instanceData->nonlinearExpressions->nl[ i]->idx;
01452                 if(foundIdx.find( index) != foundIdx.end() ){  
01453                 //if(foundIdx[ index] > 0 ){ 
01454                         //std::cout << "OLD INDEX FOUND " << index << std::endl;
01455                         //std::cout << "foundIdx[ index] " << index << std::endl;
01456                         // found an existing index
01457                         // important -- at this time m_mapExpressionTrees[ index] points to 
01458                         // the last OSExpressionTree with this index, it does not point to the 
01459                         // the just found OSExpressionTree with this index
01460                         nlNodePlus = new OSnLNodePlus();
01461                         //expTree = new OSExpressionTree(); 
01462                         expTree =  instanceData->nonlinearExpressions->nl[ i]->osExpressionTree;
01463                         // set left child to old index and right child to new one 
01464                         nlNodePlus->m_mChildren[ 0] = m_mapExpressionTrees[ index]->m_treeRoot;
01465                         nlNodePlus->m_mChildren[ 1] = instanceData->nonlinearExpressions->nl[ i]->osExpressionTree->m_treeRoot;
01466                         // we must delete the Expression tree corresponding to the old index value but not the nl nodes
01467                         instanceData->nonlinearExpressions->nl[ foundIdx[ index]  ]->m_bDeleteExpressionTree = true;
01468                         instanceData->nonlinearExpressions->nl[ foundIdx[ index]  ]->osExpressionTree->bDestroyNlNodes = false;
01469                         //point to the new expression tree
01470                         m_mapExpressionTrees[ index] = expTree;
01471                         m_mapExpressionTrees[ index]->m_treeRoot = nlNodePlus;
01472                         foundIdx[ index] = i;
01473                 }
01474                 else{  
01475                         // we have a new index
01476                         m_mapExpressionTrees[ index] = instanceData->nonlinearExpressions->nl[ i]->osExpressionTree;
01477                         m_mapExpressionTrees[ index]->m_treeRoot = instanceData->nonlinearExpressions->nl[ i]->osExpressionTree->m_treeRoot;
01478                         foundIdx[ index] = i;
01479                 }
01480                 //foundIdx[ index]++;    
01481         }
01482         // count the number of constraints and objective functions with nonlinear terms.
01483         for(pos = foundIdx.begin(); pos != foundIdx.end(); ++pos){
01484                 if(pos->first == -1) {
01485                         m_iObjectiveNumberNonlinear++; 
01486                 }
01487                 else m_iConstraintNumberNonlinear++;
01488         }
01489         foundIdx.clear();
01490         m_bProcessExpressionTrees = true;
01491         return m_mapExpressionTrees;
01492 }// getAllNonlinearExpressionTrees
01493 
01494 std::map<int, OSExpressionTree*> OSInstance::getAllNonlinearExpressionTreesMod(){
01495         if( m_bProcessExpressionTreesMod == true ) return m_mapExpressionTreesMod;
01496         m_bProcessExpressionTreesMod = true;
01497         // make sure we have the modified map available
01498         if( m_bNonLinearStructuresInitialized == false) initializeNonLinearStructures( );
01499         return m_mapExpressionTreesMod;
01500 }// getAllNonlinearExpressionTreesMod
01501 
01502 //bool OSInstance::processTimeDomain() {
01503 //      if(m_bProcessTimeDomain) return true;
01504 //      m_bProcessTimeDomain = true;
01505 //}// processTimeDomain
01506 
01507 
01508 
01509 
01510 // the set() methods
01511 
01512 bool OSInstance::setInstanceSource(string source){
01513         instanceHeader->source = source;
01514         return true;
01515 }//setInstanceSource
01516 
01517 bool OSInstance::setInstanceDescription(string description){
01518         instanceHeader->description = description;
01519         return true;
01520 }//setInstanceDescription
01521 
01522 
01523 bool OSInstance::setInstanceName(string name){
01524          instanceHeader->name = name;
01525          return true;
01526 }//setInstanceName
01527 
01528 
01529 bool OSInstance::setVariableNumber(int number){
01530         // this method assume osinstance->instanceData->variables is not null
01531         if(number <= 0) return false;
01532         //if(instanceData->variables->numberOfVariables != -1  && instanceData->variables->numberOfVariables != number){
01533         //      delete[] instanceData->variables->var;
01534         //      instanceData->variables->var = NULL;
01535         //} 
01536         if(instanceData->variables == NULL) instanceData->variables = new Variables();
01537         instanceData->variables->numberOfVariables = number;
01538         if(instanceData->variables->var == NULL){
01539                 instanceData->variables->var = new Variable*[number];
01540         }
01541         return true;
01542 }//setVariableNumber
01543 
01544 
01545 bool OSInstance::addVariable(int index, string name, double lowerBound, double upperBound, char type, double init, string initString){
01546         instanceData->variables->var[index] = new Variable();
01547         if(index < 0 || instanceData->variables->numberOfVariables <= 0 || index >= instanceData->variables->numberOfVariables) return false;
01548         instanceData->variables->var[index]->name = name;
01549         //if(lowerBound != -OSDBL_MAX && lowerBound != -OSDBL_MAX)instanceData->variables->var[index]->lb = lowerBound;
01550         instanceData->variables->var[index]->lb = lowerBound;
01551         if(upperBound != OSDBL_MAX && upperBound != OSDBL_MAX)instanceData->variables->var[index]->ub = upperBound;
01552         instanceData->variables->var[index]->type = type;
01553         if(init != OSNAN) instanceData->variables->var[index]->init = init;
01554         instanceData->variables->var[index]->initString = initString;
01555         return true;
01556 }//addVariable
01557 
01558 
01559 bool OSInstance::setVariables(int number, string *names, double *lowerBounds, 
01560         double *upperBounds, char *types, double *inits, string *initsString){
01561         if(number <= 0) return false;
01562         try{
01563                 if(instanceData->variables == NULL){
01564                         throw ErrorClass("There is no variables object");
01565                 }
01566                 if(instanceData->variables->numberOfVariables != number){
01567                         throw ErrorClass("input number of variables not equal to number in class");
01568                 }
01569                 //instanceData->variables->var = new Variable*[number];
01570                 int i;
01571                 for(i = 0; i < number; i++){
01572                         instanceData->variables->var[ i] = new Variable();
01573                 }
01574                 if(names  != NULL){
01575                         for(i = 0; i < number; i++) instanceData->variables->var[i]->name = names[i];                           
01576                 }
01577                 if(lowerBounds != NULL){
01578                         for(i = 0; i < number; i++){
01579                                 if(lowerBounds[i] != -OSDBL_MAX && lowerBounds[i] != -OSDBL_MAX)instanceData->variables->var[i]->lb = lowerBounds[i];  
01580                         }
01581                 }
01582                 if(upperBounds != NULL){
01583                         for(i = 0; i < number; i++){
01584                                 if(upperBounds[i] != OSDBL_MAX && upperBounds[i] != OSDBL_MAX)instanceData->variables->var[i]->ub = upperBounds[i]; 
01585                         }
01586                 }
01587                 if(types != NULL){
01588                         for(i = 0; i < number; i++){
01589                                 instanceData->variables->var[i]->type = types[i];
01590                                 if(types[i] != 'C' && types[i] != 'B' && types[i] != 'I' && types[i] != 'S') types[i] = 'C';
01591                         } 
01592                 }
01593                 if(inits != NULL){
01594                         for(i = 0; i < number; i++) instanceData->variables->var[i]->init = inits[i];                           
01595                 }
01596                 if(initsString != NULL){
01597                         for(i = 0; i < number; i++) instanceData->variables->var[i]->initString = initsString[i];                       
01598                 }
01599                 return true;
01600         }
01601         catch(const ErrorClass& eclass){
01602                 throw ErrorClass(  eclass.errormsg); 
01603         }
01604 }//setVariables
01605 
01606 // begin checking again with Jun Ma
01607 
01608 bool OSInstance::setObjectiveNumber(int number){
01609         if(number < 0) return false;
01610         if(instanceData->objectives == NULL) instanceData->objectives = new Objectives();
01611         if(number == 0){
01612                 instanceData->objectives->numberOfObjectives = 0;
01613                 instanceData->objectives->obj = 0;
01614                 return true;
01615         }
01616         instanceData->objectives->numberOfObjectives = number;
01617         instanceData->objectives->obj = new Objective*[number];                         
01618         return true;
01619 }//setObjectiveNumber
01620 
01621 bool OSInstance::addObjective(int index, string name, string maxOrMin, double constant, double weight, SparseVector *objectiveCoefficients){
01622         if(index >= 0 || instanceData->objectives->numberOfObjectives <= 0 || abs(index) > instanceData->objectives->numberOfObjectives) return false;
01623         int arrayIndex = abs(index) -1;
01624         if(instanceData->objectives->obj == NULL) return false;
01625         instanceData->objectives->obj[arrayIndex] = new Objective();
01626         instanceData->objectives->obj[arrayIndex]->name = name;
01627         if( (maxOrMin != "max") && (maxOrMin != "min") ) return false;
01628         else instanceData->objectives->obj[arrayIndex]->maxOrMin = maxOrMin;
01629         instanceData->objectives->obj[arrayIndex]->constant = constant;
01630         instanceData->objectives->obj[arrayIndex]->weight = weight;
01631         int n = objectiveCoefficients->number;                  
01632         instanceData->objectives->obj[arrayIndex]->numberOfObjCoef = n;
01633         if(n == 0){
01634                 instanceData->objectives->obj[arrayIndex]->coef = 0;
01635         }
01636         else{
01637                 int i = 0;
01638                 instanceData->objectives->obj[arrayIndex]->coef = new ObjCoef*[n];
01639                 for(i = 0; i < n; i++) instanceData->objectives->obj[arrayIndex]->coef[i] = new ObjCoef();
01640                 for(i = 0; i < n; i++){
01641                         instanceData->objectives->obj[arrayIndex]->coef[i]->idx = objectiveCoefficients->indexes[i];
01642                         instanceData->objectives->obj[arrayIndex]->coef[i]->value = objectiveCoefficients->values[i];   
01643                 }  
01644         }
01645         return true;
01646 }//addObjective
01647 
01648 bool OSInstance::setObjectives(int number, string *names, string *maxOrMins, double *constants, double *weights, SparseVector **objectiveCoefficients){
01649         if(number < 0) return false;
01650         try{
01651                 if(instanceData->objectives == NULL){
01652                         throw ErrorClass("there is no objectives object");              
01653                 }               
01654                 if(instanceData->objectives->numberOfObjectives != number){
01655                         throw ErrorClass("input number of objective not equal to number in class");             
01656                 }
01657                 if(number == 0) return true;
01658                 int i = 0;
01659                 for(i = 0; i < number; i++)instanceData->objectives->obj[i] = new Objective();
01660                 int j = 0;
01661                 if(names != NULL){
01662                         for(i = 0; i < number; i++) instanceData->objectives->obj[i]->name = names[i];                          
01663                 }       
01664                 if(maxOrMins != NULL){
01665                         for(i = 0; i < number; i++){
01666                                 if(maxOrMins[i] == "" || (maxOrMins[i].compare("max") != 0 && maxOrMins[i].compare("min") !=0)) return false;
01667                                 instanceData->objectives->obj[i]->maxOrMin = maxOrMins[i];                      
01668                         }
01669                 }
01670                 if(constants != NULL){
01671                         for(i = 0; i < number; i++) instanceData->objectives->obj[i]->constant = constants[i];                          
01672                 }
01673                 if(weights != NULL){
01674                         for(i = 0; i < number; i++) instanceData->objectives->obj[i]->weight = weights[i];                      
01675                 }
01676                 if(objectiveCoefficients != NULL){
01677                         for(i = 0; i < number; i++){
01678                                 int n = (&objectiveCoefficients[i] == NULL || objectiveCoefficients[i]->indexes == NULL)?0:objectiveCoefficients[i]->number;            
01679                                 instanceData->objectives->obj[i]->numberOfObjCoef = n;
01680                                 if(n == 0){
01681                                         instanceData->objectives->obj[i]->coef = NULL;
01682                                 }
01683                                 else{
01684                                         instanceData->objectives->obj[i]->coef = new ObjCoef*[n];
01685                                         for(j = 0; j < n; j++){
01686                                                 instanceData->objectives->obj[i]->coef[j] = new ObjCoef();
01687                                                 instanceData->objectives->obj[i]->coef[j]->idx  = objectiveCoefficients[i]->indexes[j];
01688                                                 instanceData->objectives->obj[i]->coef[j]->value = objectiveCoefficients[i]->values[j];                         
01689                                         }                       
01690                                 }                                                       
01691                         }
01692                 }
01693                 return true;
01694         }
01695         catch(const ErrorClass& eclass){
01696                 throw ErrorClass(  eclass.errormsg); 
01697         }
01698 }//setObjectives
01699 
01700 
01701 bool OSInstance::setConstraintNumber(int number){
01702         if(number < 0) return false;
01703         if(instanceData->constraints == NULL) instanceData->constraints = new Constraints();
01704         if(number == 0){
01705                 instanceData->constraints->numberOfConstraints = 0;
01706                 instanceData->constraints->con = 0;
01707                 return true;
01708         }
01709         instanceData->constraints->numberOfConstraints = number;
01710         if(instanceData->constraints->con == 0 ){
01711                 instanceData->constraints->con = new Constraint*[number];
01712         }
01713         return true;
01714 }//setConstraintNumber
01715 
01716 bool OSInstance::addConstraint(int index, string name, double lowerBound, double upperBound, double constant) {
01717         instanceData->constraints->con[ index] = new Constraint();
01718         if(index < 0 || instanceData->constraints->numberOfConstraints <= 0 || index >= instanceData->constraints->numberOfConstraints) return false;
01719         instanceData->constraints->con[ index]->name = name;
01720         if(lowerBound != -OSDBL_MAX && lowerBound != -OSDBL_MAX) instanceData->constraints->con[ index]->lb = lowerBound;
01721         if(upperBound != OSDBL_MAX && upperBound != OSDBL_MAX)instanceData->constraints->con[ index]->ub = upperBound;
01722         instanceData->constraints->con[ index]->constant = constant;
01723         return true;
01724 }//addConstraint
01725 
01726 
01727 bool OSInstance::setConstraints(int number, string* names, double* lowerBounds, double* upperBounds, double* constants){
01728         if(number < 0) return false;
01729         if(number == 0){
01730                 // this is done in setConstraintNumber
01731                 //instanceData->constraints = new Constraints();
01732                 //instanceData->constraints->numberOfConstraints = 0;
01733                 //instanceData->constraints->con = NULL;
01734                 return true;
01735         }
01736         try{
01737                 
01738                 if(instanceData->constraints  == NULL){
01739                         throw ErrorClass("there is no constraints object");             
01740                 }               
01741                 if(instanceData->constraints->numberOfConstraints != number){
01742                         throw ErrorClass("input number of constrasints not equal to number in class");          
01743                 }
01744                 int i = 0;
01745                 for(i = 0; i < number; i++){
01746                         instanceData->constraints->con[i] = new Constraint();
01747                 }
01748                 if(names != NULL){
01749                         for(i = 0; i < number; i++) instanceData->constraints->con[i]->name = names[i];                         
01750                 }
01751                 if(lowerBounds != NULL){
01752                         for(i = 0; i < number; i++){
01753                                 if(lowerBounds[i] != -OSDBL_MAX && lowerBounds[i] != -OSDBL_MAX)instanceData->constraints->con[i]->lb = lowerBounds[i]; 
01754                         }
01755                 }
01756                 if(upperBounds != NULL){
01757                         for(i = 0; i < number; i++){
01758                                 if(upperBounds[i] != OSDBL_MAX && upperBounds[i] != OSDBL_MAX)instanceData->constraints->con[i]->ub = upperBounds[i]; 
01759                         }
01760                 }   
01761                 if(constants != NULL){
01762                         for(i = 0; i < number; i++) instanceData->constraints->con[i]->constant = constants[i];                         
01763                 }
01764                 return true;
01765         }
01766         catch(const ErrorClass& eclass){
01767                 throw ErrorClass(  eclass.errormsg); 
01768         }
01769 }//setConstraints
01770 
01771 bool OSInstance::setLinearConstraintCoefficients(int numberOfValues, bool isColumnMajor, 
01772                 double* values, int valuesBegin, int valuesEnd, 
01773                 int* indexes, int indexesBegin, int indexesEnd,                         
01774                 int* starts, int startsBegin, int startsEnd){
01775                 if(numberOfValues < 0) return false;
01776         if(instanceData->linearConstraintCoefficients == NULL) instanceData->linearConstraintCoefficients = new LinearConstraintCoefficients() ;
01777         if(numberOfValues == 0) return true;
01778         if((values == 0 ) ||
01779            (valuesBegin < 0 || (valuesEnd - valuesBegin + 1) != numberOfValues) ||
01780            (indexes == 0) ||
01781            (indexesBegin < 0 || (indexesEnd - indexesBegin + 1) != numberOfValues) ||
01782            (starts == 0 ) ||
01783            (startsBegin < 0  || startsBegin >= startsEnd)) return false;
01784         instanceData->linearConstraintCoefficients->numberOfValues = numberOfValues;
01785         int i = 0;
01786         //starts
01787         if(instanceData->linearConstraintCoefficients->start == NULL) instanceData->linearConstraintCoefficients->start = new IntVector();
01788         if(startsBegin == 0 ){
01789                 instanceData->linearConstraintCoefficients->start->el = starts;
01790         }
01791         else{
01792                 instanceData->linearConstraintCoefficients->start->el = new int[startsEnd - startsBegin + 1];
01793                 int k = 0;
01794                 for(i = startsBegin; i <= startsEnd; i++){
01795                         instanceData->linearConstraintCoefficients->start->el[k] = starts[i];
01796                         k++;
01797                 }
01798         }       
01799         //values
01800         if(instanceData->linearConstraintCoefficients->value == NULL) instanceData->linearConstraintCoefficients->value = new DoubleVector();
01801         if(valuesBegin == 0 ){
01802                 instanceData->linearConstraintCoefficients->value->el = values;
01803         }
01804         else{
01805                 instanceData->linearConstraintCoefficients->value->el = new double[numberOfValues];
01806                 int k = 0;
01807                 for(i = valuesBegin; i <= valuesEnd; i++){
01808                         instanceData->linearConstraintCoefficients->value->el[k] = values[i];
01809                         k++;
01810                 }
01811         }
01812         //indexes
01813         if(instanceData->linearConstraintCoefficients->rowIdx == NULL) instanceData->linearConstraintCoefficients->rowIdx = new IntVector();
01814         if(instanceData->linearConstraintCoefficients->colIdx == NULL) instanceData->linearConstraintCoefficients->colIdx = new IntVector();
01815         if(isColumnMajor){
01816                 if(indexesBegin == 0 ){
01817                         instanceData->linearConstraintCoefficients->rowIdx->el = indexes;
01818                 }
01819                 else{
01820                         instanceData->linearConstraintCoefficients->rowIdx->el = new int[numberOfValues];
01821                         int k = 0;
01822                         for(i = indexesBegin; i <= indexesEnd; i++){
01823                                 instanceData->linearConstraintCoefficients->rowIdx->el[k] = indexes[i];
01824                                 k++;
01825                         }
01826                 }
01827         } 
01828         else{
01829                 if(indexesBegin == 0 ){
01830                         instanceData->linearConstraintCoefficients->colIdx->el = indexes;
01831                 }
01832                 else{
01833                         instanceData->linearConstraintCoefficients->colIdx->el = new int[numberOfValues];
01834                         int k = 0;
01835                         for(i = indexesBegin; i <= indexesEnd; i++){
01836                                 instanceData->linearConstraintCoefficients->colIdx->el[k] = indexes[i];
01837                                 k++;
01838                         }
01839                 }
01840         }
01841         return true;
01842 }//setLinearConstraintCoefficients
01843 
01844 bool OSInstance::setQuadraticTerms(int number, 
01845                 int* rowIndexes, int* varOneIndexes, int* varTwoIndexes, double* coefficients,
01846                 int begin, int end){
01847         if(number < 0) return false;
01848         if(number != (end - begin) + 1) return false;
01849         if(number == 0){
01850                 instanceData->quadraticCoefficients = 0;
01851                 return true;
01852         }
01853         if( ((end - begin + 1) != number) ||
01854            (rowIndexes == 0) ||                    
01855            (varOneIndexes == 0) ||
01856            (varTwoIndexes == 0) ||
01857            (coefficients == 0) ) return false;
01858         instanceData->quadraticCoefficients = new QuadraticCoefficients();
01859         instanceData->quadraticCoefficients->numberOfQuadraticTerms = number;
01860         int i = 0;
01861         instanceData->quadraticCoefficients->qTerm = new QuadraticTerm*[number];
01862         for(i = 0; i < number; i++) instanceData->quadraticCoefficients->qTerm[i] = new QuadraticTerm();
01863         int k = 0;
01864         for(i = begin; i <= end; i++){
01865                         instanceData->quadraticCoefficients->qTerm[k]->idx = rowIndexes[i];
01866                         instanceData->quadraticCoefficients->qTerm[k]->idxOne = varOneIndexes[i];
01867                         instanceData->quadraticCoefficients->qTerm[k]->idxTwo = varTwoIndexes[i];
01868                         instanceData->quadraticCoefficients->qTerm[k]->coef = coefficients[i];
01869                         k++;
01870         }
01871         return true;
01872 }//setQuadraticTerms
01873 
01874 bool OSInstance::setQuadraticTermsInNonlinearExpressions(int numQPTerms, int* rowIndexes, int* varOneIndexes, int* varTwoIndexes, double* coefficients){
01875                 instanceData->nonlinearExpressions->numberOfNonlinearExpressions = numQPTerms;
01876                 instanceData->nonlinearExpressions->nl = new Nl*[ numQPTerms ];
01877                 // define the vectors
01878                 OSnLNode *nlNodePoint;
01879                 OSnLNodeVariable *nlNodeVariablePoint;
01880                 std::vector<OSnLNode*> nlNodeVec;
01881                 //
01882                 //
01883                 int i;
01884                 for(i = 0; i < numQPTerms; i++){
01885                         instanceData->nonlinearExpressions->nl[ i] = new Nl();
01886                         instanceData->nonlinearExpressions->nl[ i]->idx = rowIndexes[ i];
01887                         instanceData->nonlinearExpressions->nl[ i]->osExpressionTree = new OSExpressionTree();
01888                 // create a variable nl node for x0
01889                 nlNodeVariablePoint = new OSnLNodeVariable();
01890                 nlNodeVariablePoint->idx = varOneIndexes[ i];
01891                 // give this variable the coefficient
01892                 nlNodeVariablePoint->coef = coefficients[ i];
01893                 nlNodeVec.push_back( nlNodeVariablePoint);
01894                 // create the nl node for x1
01895                 nlNodeVariablePoint = new OSnLNodeVariable();
01896                 nlNodeVariablePoint->idx = varTwoIndexes[ i];
01897                 nlNodeVec.push_back( nlNodeVariablePoint);
01898                 // create the nl node for *
01899                 nlNodePoint = new OSnLNodeTimes();
01900                 nlNodeVec.push_back( nlNodePoint);
01901                 // the vectors are in postfix format
01902                 // now the expression tree
01903                 instanceData->nonlinearExpressions->nl[ i]->osExpressionTree->m_treeRoot =
01904                         nlNodeVec[ 0]->createExpressionTreeFromPostfix( nlNodeVec);
01905                 nlNodeVec.clear();
01906                 }
01907         return true;
01908 }//setQuadraticTermsInNonlinearExpressions
01909 
01910 bool OSInstance::initializeNonLinearStructures( ){
01911         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
01912         if( m_bNonLinearStructuresInitialized == true) return true;
01913         if( m_bProcessVariables == false) processVariables();
01914         if( m_bProcessObjectives == false) processObjectives();
01915         if( m_bProcessConstraints == false) processConstraints();
01916         m_iVariableNumber = instanceData->variables->numberOfVariables;
01917         m_iConstraintNumber = instanceData->constraints->numberOfConstraints;
01918         m_iObjectiveNumber = instanceData->objectives->numberOfObjectives;
01919         // get all of the expression trees
01920         if( m_bProcessExpressionTrees == false) getAllNonlinearExpressionTrees();
01921         // before proceeding get a copy of the map of the Expression Trees
01922         if( m_bDuplicateExpressionTreesMap == false) duplicateExpressionTreesMap();
01923         // now create all of the variable maps for each expression tree
01924         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){  
01925                 (posMapExpTree->second)->getVariableIndiciesMap() ;
01926         }
01927         // add the quadratic terms if necessary
01928         if(getNumberOfQuadraticTerms() > 0) addQTermsToExressionTree();
01929         // now get the map of all nonlinear variables
01930         getAllNonlinearVariablesIndexMap( );
01931         getDenseObjectiveCoefficients();
01932         m_mdConstraintFunctionValues = new double[ this->instanceData->constraints->numberOfConstraints];
01933         m_mdObjectiveFunctionValues = new double[ this->instanceData->objectives->numberOfObjectives];
01934         //m_mdObjGradient = new double[ this->instanceData->variables->numberOfVariables];
01935         m_bNonLinearStructuresInitialized = true;
01936         return true;
01937 }
01938 
01939 SparseJacobianMatrix *OSInstance::getJacobianSparsityPattern( ){
01940         // if already called return the sparse Jacobian
01941         // it is important that this method NOT get called twice -- if
01942         // there are linear terms in <linearConstraintCoefficients> that
01943         // also appear in <nonlinearExpressions> then they will keep getting added
01944         // to the modified expession tree with each call to this method
01945         if( m_bSparseJacobianCalculated == true) return m_sparseJacMatrix;
01946         //std::cout << "INSIDE GET JACOBIAN SPARSITY PATTERN" << std::endl;
01947         // determine if we are in column or row major
01948         getLinearConstraintCoefficientMajor();
01949         // make sure the data structures have been inialized
01950         if( m_bNonLinearStructuresInitialized == false) initializeNonLinearStructures( );
01951         try{
01952                 if( m_bColumnMajor == true){
01953                          if( getSparseJacobianFromColumnMajor( ) == false) throw ErrorClass("An error occurred in getSpareJacobianFromColumnMajor");
01954                 }
01955                 else {
01956                         if( getSparseJacobianFromRowMajor( ) == false) throw ErrorClass("An error occurred in getSpareJacobianFromRowMajor");
01957                 }
01958         }
01959         catch(const ErrorClass& eclass){
01960                 throw ErrorClass(  eclass.errormsg); 
01961         }
01962         // now fill in the arrays of the sparseJacMatrix
01963         m_sparseJacMatrix = new SparseJacobianMatrix();
01964         // we point to memory already created so don't 
01965         // destroy during garbage collection
01966         m_sparseJacMatrix->bDeleteArrays = false;
01967         m_sparseJacMatrix->valueSize =  m_iJacValueSize;
01968         m_sparseJacMatrix->starts = m_miJacStart;
01969         m_sparseJacMatrix->conVals = m_miJacNumConTerms;
01970         m_sparseJacMatrix->indexes = m_miJacIndex;
01971         m_sparseJacMatrix->values = m_mdJacValue;
01972         m_bSparseJacobianCalculated = true;
01973         return m_sparseJacMatrix;
01974 }//getJacobianSparsityPatter
01975 
01976 bool OSInstance::addQTermsToExressionTree(){
01977         int i, k, idx;
01978         // get the number of qTerms
01979         int numQTerms = instanceData->quadraticCoefficients->numberOfQuadraticTerms;    
01980         if(numQTerms <= 0 || m_bQTermsAdded == true) return true;
01981         OSnLNodeVariable* nlNodeVariableOne;
01982         OSnLNodeVariable* nlNodeVariableTwo;
01983         OSnLNodeTimes* nlNodeTimes;
01984         OSnLNodePlus* nlNodePlus;
01985         OSExpressionTree* expTree;
01986         getQuadraticTerms();    
01987         std::cout << "PROCESSING QUADRATIC TERMS" << std::endl;
01988         for(i = 0; i < numQTerms; i++){
01989                 idx = m_quadraticTerms->rowIndexes[ i];
01990                 std::cout << "PROCESSING QTERM  = "  << i <<std::endl;
01991                 // see if row idx is in the expression tree
01992                 if( m_mapExpressionTreesMod.find( idx) != m_mapExpressionTreesMod.end() ) {
01993                         // row idx is in the expression tree
01994                         // add the qTerm in row idx  to the expression tree     
01995                         // define two new OSnLVariable nodes, an OSnLnodeTimes, and OSnLnodePlus 
01996                         nlNodeVariableOne = new OSnLNodeVariable();
01997                         nlNodeVariableOne->idx = m_quadraticTerms->varOneIndexes[ i];
01998                         // see if the variable indexed by nlNodeVariableOne->idx is in the expression tree for row idx
01999                         // if not, add to mapVarIdx
02000                         expTree = m_mapExpressionTreesMod[ idx];
02001                         if(  expTree->m_bIndexMapGenerated == false) expTree->getVariableIndiciesMap();
02002                         if( (*expTree->mapVarIdx).find( nlNodeVariableOne->idx) == (*expTree->mapVarIdx).end()  ){
02003                                 // add to map
02004                                 k = (*expTree->mapVarIdx).size();
02005                                 (*expTree->mapVarIdx)[ nlNodeVariableOne->idx] =  k + 1;
02006                                 std::cout << "ADDED THE FOLLOWING VAIRABLE TO THE MAP" <<  nlNodeVariableOne->idx << std::endl;
02007                         }
02008                         nlNodeVariableOne->coef = m_quadraticTerms->coefficients[ i];
02009                         nlNodeVariableTwo = new OSnLNodeVariable();
02010                         nlNodeVariableTwo->idx = m_quadraticTerms->varTwoIndexes[ i];
02011                         // see if the variable indexed by nlNodeVariableTwo->idx is in the expression tree for row idx
02012                         // if not, add to mapVarIdx
02013                         if( (*expTree->mapVarIdx).find( nlNodeVariableTwo->idx) == (*expTree->mapVarIdx).end()  ){
02014                                 // add to map
02015                                 k = (*expTree->mapVarIdx).size();
02016                                 (*expTree->mapVarIdx)[ nlNodeVariableTwo->idx] =  k + 1;
02017                                 std::cout << "ADDED THE FOLLOWING VAIRABLE TO THE MAP" <<  nlNodeVariableTwo->idx << std::endl;
02018                         }
02019                         nlNodeVariableTwo->coef = 1.;
02020                         // now multiply the two new variable nodes together
02021                         nlNodeTimes = new OSnLNodeTimes();
02022                         nlNodeTimes->m_mChildren[ 0] = nlNodeVariableOne;
02023                         nlNodeTimes->m_mChildren[ 1] = nlNodeVariableTwo;               
02024                         // now add the result to the expression tree
02025                         nlNodePlus = new OSnLNodePlus();
02026                         nlNodePlus->m_mChildren[ 0] = expTree->m_treeRoot;
02027                         nlNodePlus->m_mChildren[ 1] = nlNodeTimes;
02028                         //expTree = new OSExpressionTree();
02029                         expTree->m_treeRoot = nlNodePlus ;
02030                         // get rid of old variable map
02031                         if(expTree->m_bIndexMapGenerated == true){
02032                                 delete expTree->mapVarIdx;
02033                                 expTree->mapVarIdx = NULL;
02034                                 expTree->m_bIndexMapGenerated = false;
02035                         }       
02036                         //expTree->m_bIndexMapGenerated = false;
02037                         //m_mapExpressionTreesMod[ idx ]  = expTree;    
02038                         //expTree->mapVarIdx = m_mapExpressionTreesMod[ idx]->mapVarIdx;
02039                 }
02040                 else{ 
02041                         // create the quadratic expression to add to the expression tree
02042                         nlNodeVariableOne = new OSnLNodeVariable();
02043                         nlNodeVariableOne->idx = m_quadraticTerms->varOneIndexes[ i];
02044                         nlNodeVariableOne->coef = m_quadraticTerms->coefficients[ i];
02045                         nlNodeVariableTwo = new OSnLNodeVariable();
02046                         nlNodeVariableTwo->idx = m_quadraticTerms->varTwoIndexes[ i];
02047                         nlNodeVariableTwo->coef = 1.;
02048                         // now multiply the two new variable nodes together
02049                         nlNodeTimes = new OSnLNodeTimes();
02050                         nlNodeTimes->m_mChildren[ 0] = nlNodeVariableOne;
02051                         nlNodeTimes->m_mChildren[ 1] = nlNodeVariableTwo;
02052                         // create a new expression tree corresponding to row idx.
02053                         expTree = new OSExpressionTree();                                               
02054                         expTree->m_treeRoot = nlNodeTimes ;
02055                         expTree->mapVarIdx = expTree->getVariableIndiciesMap();         
02056                         m_mapExpressionTreesMod[ idx ]  = expTree;
02057                         if(idx < 0){
02058                                 m_iObjectiveNumberNonlinear++;
02059                                 m_bProcessExpressionTrees = true;
02060                         }
02061                         else{
02062                                 m_iConstraintNumberNonlinear++;
02063                                 m_bProcessExpressionTrees = true;
02064                         }
02065                         std::cout << "NUMBER OF EXPRESSION TREES = "  << m_mapExpressionTreesMod.size() <<std::endl;
02066                         std::cout << "NUMBER OF NONLINEAR OBJECTIVES = "  << getNumberOfNonlinearObjectives() <<std::endl;
02067                 } 
02068                 // if there were no nonlinear terms make this the expression tree
02069                 if(m_iNonlinearExpressionNumber <= 0) m_mapExpressionTrees = m_mapExpressionTreesMod;
02070                 m_bQTermsAdded =true;
02071         }
02072         return true;
02073 }
02074 
02075 double OSInstance::calculateFunctionValue(int idx, double *x, bool new_x){
02076         try{
02077 
02078                 int i, j;
02079                 double dvalue = 0;
02080                 if( m_binitForAlgDiff == false) initForAlgDiff();
02081                 if( m_bSparseJacobianCalculated == false) getJacobianSparsityPattern();
02082                 if(idx >= 0){ // we have a constraint
02083                         // make sure the index idx is valid
02084                         if( getConstraintNumber() <= idx  ) throw 
02085                         ErrorClass("constraint index not valid in OSInstance::calculateFunctionValue");
02086                         if( new_x == false) return *(m_mdConstraintFunctionValues + idx);
02087                         // get the nonlinear part
02088                         if( m_mapExpressionTreesMod.find( idx) != m_mapExpressionTreesMod.end() ){
02089                                 dvalue = m_mapExpressionTreesMod[ idx]->calculateFunction( x,  new_x);
02090                                 //dvalue = vdFunVals[ idx + 1];
02091                         }
02092                         // now the linear part
02093                         // be careful, loop over only the constant terms in sparseJacMatrix
02094                         i = m_sparseJacMatrix->starts[ idx];
02095                         j = m_sparseJacMatrix->starts[ idx + 1 ];
02096                         //while ( i <  j &&  (i - m_sparseJacMatrix->starts[ idx])  < m_sparseJacMatrix->conVals[ idx] ){
02097                         while ( (i - m_sparseJacMatrix->starts[ idx])  < m_sparseJacMatrix->conVals[ idx] ){
02098                                 //std::cout << "m_sparseJacMatrix->values[ i] " << m_sparseJacMatrix->values[ i] << std::endl;
02099                                 //std::cout << "m_sparseJacMatrix->indexes[ i] " << m_sparseJacMatrix->indexes[ i] << std::endl;
02100                                 dvalue += m_sparseJacMatrix->values[ i]*x[ m_sparseJacMatrix->indexes[ i] ];
02101                                 i++;
02102                         }       
02103                         // add in the constraint function constant
02104                         dvalue += m_mdConstraintConstants[ idx ];
02105                         return dvalue;
02106                 }
02107                 else{ // we have an objective function
02108                         // make sure the index idx is valid
02109                         if( getObjectiveNumber() <= ( abs( idx) - 1) ) throw 
02110                         ErrorClass("objective function index not valid in OSInstance::calculateFunctionValue");
02111                         if( new_x == false) return *(m_mdObjectiveFunctionValues + ( abs( idx) - 1));
02112                         // get the nonlinear part
02113                         if( m_mapExpressionTreesMod.find( idx) != m_mapExpressionTreesMod.end() ){
02114                                 dvalue = m_mapExpressionTreesMod[ idx]->calculateFunction( x,  new_x);
02115                         }
02116                         // get linear part
02117                         SparseVector **objCoef = getObjectiveCoefficients();
02118                         SparseVector *obj = objCoef[ abs( idx) - 1];
02119                         for(i = 0; i < obj->number; i++){
02120                                 dvalue += x[ obj->indexes[i]]*(obj->values[ i]);
02121                         }
02122                         // add in the objective function constant
02123                         dvalue += m_mdObjectiveConstants[ abs( idx) - 1 ];
02124                         // get the coefficients for objective function idx
02125                         *(m_mdObjectiveFunctionValues + ( abs( idx) - 1)) = dvalue;
02126                         return *(m_mdObjectiveFunctionValues + ( abs( idx) - 1));
02127                 }
02128         }
02129         catch(const ErrorClass& eclass){
02130                 throw ErrorClass( eclass.errormsg);
02131         } 
02132 }//calculateFunctionValue
02133 
02134 
02135 double *OSInstance::calculateAllConstraintFunctionValues( double* x, double *objLambda, double *conLambda,
02136         bool new_x, int highestOrder){  
02137         try{
02138                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) 
02139                         getIterateResults(x, objLambda, conLambda, new_x,  highestOrder);
02140         }
02141         catch(const ErrorClass& eclass){
02142                 throw ErrorClass( eclass.errormsg);
02143         } 
02144         return m_mdConstraintFunctionValues;
02145 }//calculateAllConstraintFunctionValues
02146 
02147 
02148 double *OSInstance::calculateAllConstraintFunctionValues(double* x, bool new_x){
02149         try{
02150                 m_iHighestOrderEvaluated = -1;
02151                 if( new_x == false) return m_mdConstraintFunctionValues;
02152                 int idx, numConstraints;
02153                 numConstraints = getConstraintNumber();
02154                 // loop over all constraints
02155                 for(idx = 0; idx < numConstraints; idx++){
02156                         m_mdConstraintFunctionValues[ idx]  = calculateFunctionValue(idx, x, new_x);    
02157                 }
02158                 
02159         }
02160         catch(const ErrorClass& eclass){
02161                 throw ErrorClass( eclass.errormsg);
02162         } 
02163         return m_mdConstraintFunctionValues;    
02164 }//end calculateAllConstraintFunctionValues
02165 
02166 
02167 double *OSInstance::calculateAllObjectiveFunctionValues( double* x, double *objLambda, double *conLambda,
02168         bool new_x, int highestOrder){  
02169         try{
02170                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) 
02171                         getIterateResults(x, objLambda, conLambda, new_x,  highestOrder);
02172         }
02173         catch(const ErrorClass& eclass){
02174                 throw ErrorClass( eclass.errormsg);
02175         } 
02176         return m_mdObjectiveFunctionValues;
02177 }//calculateAllConstraintFunctionValues
02178 
02179 
02180 double *OSInstance::calculateAllObjectiveFunctionValues( double* x, bool new_x){        
02181         try{
02182                 m_iHighestOrderEvaluated = -1;
02183                 if( new_x == false) return m_mdObjectiveFunctionValues;
02184                 int idx, numObjectives;
02185                 numObjectives = getObjectiveNumber();
02186                 // loop over all constraints
02187                 for(idx = 0; idx < numObjectives; idx++){
02188                         m_mdObjectiveFunctionValues[ idx]  = calculateFunctionValue(-idx -1, x, new_x); 
02189                 }
02190         }
02191         catch(const ErrorClass& eclass){
02192                 throw ErrorClass( eclass.errormsg);
02193         } 
02194         return m_mdObjectiveFunctionValues;
02195 }//calculateAllObjectiveFunctionValues
02196 
02197 
02198 SparseJacobianMatrix *OSInstance::calculateAllConstraintFunctionGradients(double* x, double *objLambda, double *conLambda,
02199                 bool new_x, int highestOrder){
02200         try{
02201                 if(highestOrder < 1 ) throw ErrorClass("When calling calculateAllConstraintFunctionGradients highestOrder should be 1 or 2");
02202                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) 
02203                         getIterateResults(x, objLambda, conLambda,  new_x,  highestOrder);
02204         }//end try
02205         catch(const ErrorClass& eclass){
02206                 throw ErrorClass( eclass.errormsg);
02207         } 
02208         return m_sparseJacMatrix;
02209 }//calculateAllConstraintFunctionGradients      
02210 
02211 
02212 
02213 SparseVector *OSInstance::calculateConstraintFunctionGradient(double* x, double *objLambda, double *conLambda,
02214                 int idx, bool new_x, int highestOrder){
02215         try{
02216                 if(highestOrder < 1 ) throw ErrorClass("When calling calculateConstraintFunctionGradient highestOrder should be 1 or 2");
02217                 if(idx > instanceData->variables->numberOfVariables ) 
02218                         throw ErrorClass("invalid index passed to calculateConstraintFunctionGrad");
02219                 SparseVector *sp;
02220                 sp = new SparseVector();
02221                 sp->bDeleteArrays = true;
02222                 int i;
02223                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) 
02224                         getIterateResults(x, objLambda, conLambda,  new_x,  highestOrder);
02225                 sp->number = m_miJacStart[ idx + 1] - m_miJacStart[ idx];
02226                 sp->values = new double[ sp->number];
02227                 sp->indexes = new int[ sp->number];
02228                 for(i = 0; i < sp->number; i++){
02229                         sp->values[ i] = m_mdJacValue[ m_miJacStart[ idx] +  i];
02230                         sp->indexes[ i] = m_miJacIndex[ m_miJacStart[ idx] +  i];
02231                 }
02232                 return sp;
02233         }
02234         catch(const ErrorClass& eclass){
02235                 throw ErrorClass( eclass.errormsg);
02236         } 
02237 }//calculateConstraintFunctionGradient
02238 
02239 
02240 SparseVector *OSInstance::calculateConstraintFunctionGradient(double* x, int idx, bool new_x){
02241         try{
02242                 if(idx > instanceData->variables->numberOfVariables ) 
02243                         throw ErrorClass("invalid index passed to calculateConstraintFunctionGrad");
02244                 SparseVector *sp;
02245                 sp = new SparseVector();
02246                 sp->bDeleteArrays = true;
02247                 int i;
02248                 if( new_x == true || (1 > m_iHighestOrderEvaluated)  ) 
02249                         getIterateResults(x, NULL, NULL,  new_x,  1);
02250                 sp->number = m_miJacStart[ idx + 1] - m_miJacStart[ idx];
02251                 sp->values = new double[ sp->number];
02252                 sp->indexes = new int[ sp->number];
02253                 for(i = 0; i < sp->number; i++){
02254                         sp->values[ i] = m_mdJacValue[ m_miJacStart[ idx] +  i];
02255                         sp->indexes[ i] = m_miJacIndex[ m_miJacStart[ idx] + i];
02256                 }
02257                 return sp;
02258         }
02259         catch(const ErrorClass& eclass){
02260                 throw ErrorClass( eclass.errormsg);
02261         } 
02262 }//calculateConstraintFunctionGradient
02263 
02264 
02265 
02266 double **OSInstance::calculateAllObjectiveFunctionGradients(double* x, double *objLambda, double *conLambda,
02267                 bool new_x, int highestOrder){
02268         try{
02269                 if(highestOrder < 1 ) throw ErrorClass("When calling calculateAllObjectiveFunctionGradients highestOrder should be 1 or 2");
02270                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) {
02271                         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02272                         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02273                                 if(posMapExpTree->first < 0){ // this nonlinear expression indexes an objective function
02274                                         m_mmdObjGradient[ abs( posMapExpTree->first) - 1 ] = calculateObjectiveFunctionGradient(x, objLambda, conLambda, 
02275                                                         posMapExpTree->first, new_x, highestOrder);
02276                                 }       
02277                         }
02278                 }
02279         }
02280         catch(const ErrorClass& eclass){
02281                 throw ErrorClass( eclass.errormsg);
02282         } 
02283         return m_mmdObjGradient;
02284 }// calculateAllObjectiveFunctionGradients
02285 
02286 double *OSInstance::calculateObjectiveFunctionGradient(double* x, double *objLambda, double *conLambda,
02287                 int objIdx, bool new_x, int highestOrder){
02288         /* if we are just doing an objective function gradient we should do a zero order forward sweep
02289          * and a reverse first order reverse sweep
02290          */
02291         try{
02292                 if(highestOrder < 1 ) throw ErrorClass("When calling calculateObjectiveFunctionGradient highestOrder should be 1 or 2");        
02293                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) {
02294                         int domainIdx = 0;      
02295                         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02296                         std::map<int, int>::iterator posVarIndexMap;
02297                         unsigned int i;
02298                         int iHighestOrderEvaluatedStore;        
02299                         iHighestOrderEvaluatedStore = m_iHighestOrderEvaluated;
02300                         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02301                                 //kipp: modify for more than one obj
02302                                 if(posMapExpTree->first == objIdx){
02303                                         if( new_x == true ){
02304                                                 if( m_vdX.size() > 0) m_vdX.clear();
02305                                                 for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
02306                                                         m_vdX.push_back( x[ posVarIndexMap->first]) ;
02307                                                 }
02308                                                 if( (m_bCppADFunIsCreated == false || m_bCppADMustReTape == true )  && (m_mapExpressionTreesMod.size() > 0) ) {
02309                                                         createCppADFun( m_vdX);
02310                                                 }                                               
02311                                         }       
02312                                         
02313                                         if(( new_x == true ) || (m_iHighestOrderEvaluated < 0) )this->forwardAD(0, m_vdX);
02314                                         
02315                                         if(( new_x == true ) || (m_iHighestOrderEvaluated < 1) ) {
02316                                                 m_vdRangeUnitVec[ domainIdx] = 1.;
02317                                                 m_vdYjacval = this->reverseAD(1, m_vdRangeUnitVec);
02318                                                 for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
02319                                                         m_mmdObjGradient[ abs( objIdx) - 1 ][ m_miNonLinearVarsReverseMap[ i]] = m_vdYjacval[ i] + 
02320                                                                 m_mmdDenseObjectiveCoefficients[  abs( objIdx) - 1][ m_miNonLinearVarsReverseMap[ i]];
02321                                                 }
02322                                         }
02323                                         m_iHighestOrderEvaluated = iHighestOrderEvaluatedStore;
02324                                         m_vdRangeUnitVec[ domainIdx] = 0.;
02325                                         // exit the loop
02326                                         break;
02327                                 }
02328                                 domainIdx++;
02329                         }
02330                 }                       
02331         }
02332         catch(const ErrorClass& eclass){
02333                 throw ErrorClass( eclass.errormsg);
02334         } 
02335         return m_mmdObjGradient[abs( objIdx) - 1];
02336 }//calculateObjectiveFunctionGradient
02337 
02338 
02339 double *OSInstance::calculateObjectiveFunctionGradient(double* x, int objIdx, bool new_x){
02340         try{
02341                 int domainIdx = 0;      
02342                 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02343                 std::map<int, int>::iterator posVarIndexMap;
02344                 unsigned int i;
02345                 int iHighestOrderEvaluatedStore;        
02346                 iHighestOrderEvaluatedStore = m_iHighestOrderEvaluated;
02347                 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02348                         if(posMapExpTree->first == objIdx){
02349                                 if( new_x == true ){
02350                                         if( m_vdX.size() > 0) m_vdX.clear();
02351                                         for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
02352                                                 m_vdX.push_back( x[ posVarIndexMap->first]) ;
02353                                         }
02354                                         if( (m_bCppADFunIsCreated == false || m_bCppADMustReTape == true )  && (m_mapExpressionTreesMod.size() > 0) ) {
02355                                                 createCppADFun( m_vdX);
02356                                         }                                               
02357                                 }       
02358                                 
02359                                 if(( new_x == true ) || (m_iHighestOrderEvaluated < 0) )this->forwardAD(0, m_vdX);
02360                                 
02361                                 if(( new_x == true ) || (m_iHighestOrderEvaluated < 1) ) {
02362                                         m_vdRangeUnitVec[ domainIdx] = 1.;
02363                                         m_vdYjacval = this->reverseAD(1, m_vdRangeUnitVec);
02364                                         for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
02365                                                 m_mmdObjGradient[ abs( objIdx) - 1 ][ m_miNonLinearVarsReverseMap[ i]] = m_vdYjacval[ i] + 
02366                                                         m_mmdDenseObjectiveCoefficients[  abs( objIdx) - 1][ m_miNonLinearVarsReverseMap[ i]];
02367                                         }
02368                                 }
02369                                 m_iHighestOrderEvaluated = iHighestOrderEvaluatedStore;
02370                                 m_vdRangeUnitVec[ domainIdx] = 0.;
02371                                 // exit the loop
02372                                 break;
02373                         }
02374                         domainIdx++;
02375                 }
02376         }       
02377         catch(const ErrorClass& eclass){
02378                 throw ErrorClass( eclass.errormsg);
02379         } 
02380         return m_mmdObjGradient[abs( objIdx) - 1];
02381 }//calculateObjectiveFunctionGradient   
02382 
02383 
02384 SparseHessianMatrix *OSInstance::calculateLagrangianHessian( double* x, double *objLambda, double *conLambda,
02385         bool new_x, int highestOrder){
02386         try{
02387                 if(highestOrder != 2 ) throw ErrorClass("When calling calculateLagrangianHessian highestOrder should be 2");
02388                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) {
02389                         //std::cout  << "CALL getIterateResults() FROM calculateLagrangianHessain" << std::endl;
02390                         getIterateResults(x, objLambda, conLambda,  new_x,  highestOrder);
02391                 }
02392         }
02393         catch(const ErrorClass& eclass){
02394                 throw ErrorClass( eclass.errormsg);
02395         } 
02396         return m_LagrangianSparseHessian;
02397 }//calculateLagrangianHessian
02398 
02399 SparseHessianMatrix *OSInstance::calculateHessian(double* x, int idx, bool new_x){
02400         try{
02401                 if(idx > instanceData->variables->numberOfVariables ) 
02402                         throw ErrorClass("invalid index passed to calculateHessian");
02403                 double *objLambda = new double[ getObjectiveNumber() ];
02404                 double *conLambda = new double[ getConstraintNumber() ];
02405                 //std::cout << "NUMBER OF OBJECTIVES = " << getObjectiveNumber() << std::endl;
02406                 //std::cout << "NUMBER OF CONSTRAINTS = " << getConstraintNumber() << std::endl;
02407                 int i;
02408                 // initialize all to zero
02409                 for(i = 0; i < getObjectiveNumber(); i++){
02410                         objLambda[ i] = 0.0;
02411                 }
02412                 for(i = 0; i < getConstraintNumber(); i++){
02413                         conLambda[ i] = 0.0;
02414                 }
02415                 // see if we have the index of an objective function or a constraint
02416                 // and more corresponding component 1.0
02417                 if(idx < 0){
02418                         objLambda[ abs(idx) - 1] = 1.0;
02419                 }
02420                 else{
02421                         conLambda[ idx] = 1.0;
02422                 }
02423                 if( new_x == true || (2 > m_iHighestOrderEvaluated)  ) {
02424                         getIterateResults(x, objLambda, conLambda, new_x,  2);
02425                         std::cout  << "CALL getIterateResults() FROM calculateHessian" << std::endl;
02426                 }
02427                 delete[] objLambda;
02428                 delete[] conLambda;
02429         }
02430         catch(const ErrorClass& eclass){
02431                 throw ErrorClass( eclass.errormsg);
02432         } 
02433         return m_LagrangianSparseHessian;
02434 }//calculateHessian                     
02435 
02436 
02437 
02438 bool OSInstance::getSparseJacobianFromColumnMajor( ){
02439         // we assume column major matrix
02440         if( m_bColumnMajor == false) return false;
02441         int iNumRowStarts = getConstraintNumber() + 1;  
02442         int i,j, iTemp;
02443         int iNumVariableStarts = getVariableNumber() ;
02444         int *start = NULL;
02445         int *index = NULL;
02446         double *value = NULL;
02447         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02448                 start = this->instanceData->linearConstraintCoefficients->start->el;
02449                 index = this->instanceData->linearConstraintCoefficients->rowIdx->el;
02450                 value = this->instanceData->linearConstraintCoefficients->value->el;
02451         }
02452         m_miJacStart = new int[ iNumRowStarts];
02453         m_miJacNumConTerms = new int[ getConstraintNumber()];
02454         OSnLNodePlus *nlNodePlus;
02455         OSnLNodeVariable *nlNodeVariable;
02456         OSExpressionTree *expTree = NULL;
02457         // now initialize starts and variable index maps 
02458         for ( i = 0; i < iNumRowStarts; i++){                   
02459                 m_miJacStart [ i ] = 0;
02460                 // map the variables  in the nonlinear rows
02461                 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
02462                         // the following is equivalent to  m_treeRoot->getVariableIndexMap( i);
02463                         m_mapExpressionTreesMod[ i]->getVariableIndiciesMap();
02464                         
02465                 }
02466         }
02467         // only execute the following code if there are linear constraint coefficients
02468         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02469                 // i is indexing columns (variables) and j is indexing row numbers 
02470                 for (i = 0; i < iNumVariableStarts; i++){       
02471                         for (j = start[i]; j < start[ i + 1 ]; j++){
02472                                 // index[ j] is a row index, we have just found an occurance of row index[j]
02473                                 // therefore we increase by 1 (or push back) the start of the row indexed by index[j] + 1, 
02474                                 // i.e. the start of the next row
02475                                 // check to see if variable i is linear/constant in the row index[ j] 
02476                                 // if so, increment m_miJacStart[ index[j] + 1]
02477                                 //
02478                                 if( (m_mapExpressionTreesMod.find( index[ j]) != m_mapExpressionTreesMod.end() ) &&
02479                                         ( (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).find( i) != (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).end()) ){
02480                                         // variable i is appears in the expression tree for row index[ j]
02481                                         // add the coefficient corresponding to variable i in row index[ j] to the expression tree      
02482                                         // define a new OSnLVariable and OSnLnodePlus 
02483                                         expTree = m_mapExpressionTreesMod[ index[j]  ];
02484                                         nlNodeVariable = new OSnLNodeVariable();
02485                                         nlNodeVariable->coef = value[ j];
02486                                         nlNodeVariable->idx = i;
02487                                         nlNodePlus = new OSnLNodePlus();
02488                                         nlNodePlus->m_mChildren[ 0] = m_mapExpressionTreesMod[ index[ j] ]->m_treeRoot;
02489                                         nlNodePlus->m_mChildren[ 1] = nlNodeVariable;
02490                                         //expTree = new OSExpressionTree();
02491                                         expTree->m_treeRoot = nlNodePlus ;      
02492                                         //expTree->mapVarIdx = m_mapExpressionTreesMod[ index[ j]]->mapVarIdx;
02493                                         //m_mapExpressionTreesMod[ index[ j] ]  = expTree;
02494                                         //std::cout << m_mapExpressionTreesMod[ index[ j] ]->m_treeRoot->getNonlinearExpressionInXML() << std::endl;    
02495                                         //std::cout << m_mapExpressionTrees[ index[ j] ]->m_treeRoot->getNonlinearExpressionInXML() << std::endl;
02496                                 }
02497                                 else{ 
02498                                         m_miJacStart[ index[j] + 1] ++;
02499                                 }                               
02500                         }
02501                 }
02502         }
02503         // at this point, m_miJacStart[ i] holds the number of columns with a linear/constant nonzero in row i - 1
02504         // we are not done with the start indicies, if we are here, and we
02505         // knew the correct starting point of row i -1, the correct starting point
02506         // for row i is m_miJacStart[i] + m_miJacStart [i - 1]
02507         m_miJacStart[0] = 0;
02508         for (i = 1; i < iNumRowStarts; i++ ){
02509                 m_miJacNumConTerms[ i - 1] = m_miJacStart[i];
02510                 if( m_mapExpressionTreesMod.find( i - 1) != m_mapExpressionTreesMod.end() ){
02511                         m_miJacStart[i] += (m_miJacStart[i - 1] + (*m_mapExpressionTreesMod[ i - 1]->mapVarIdx).size() );
02512                 }
02513                 else{
02514                         m_miJacStart[i] += m_miJacStart[i - 1];
02515                 }       
02516         }
02517         // dimension miIndex and mdValue here
02518         m_iJacValueSize =       m_miJacStart[ iNumRowStarts - 1];
02519         m_miJacIndex = new int[  m_iJacValueSize];
02520         m_mdJacValue = new double[ m_iJacValueSize ];
02521         // now get the values of the constant terms if there are any
02522         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02523                 // loop over variables  
02524                 for (i = 0; i < iNumVariableStarts; i++){
02525                         // get row indices and values of the A matrix
02526                         // kipp -- should we have a check to see if start[i+1] > start[i]
02527                         for (j = start[i]; j < start[ i + 1 ]; j++){
02528                                 // store this variable index in every row where the variable appears
02529                                 // however, don't store this as constant term if it appears in mapVarIdx
02530                                 if( (m_mapExpressionTreesMod.find( index[ j]) == m_mapExpressionTreesMod.end() ) ||
02531                                         ( (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).find( i) == (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).end()) ){
02532                                         iTemp = m_miJacStart[ index[j]];
02533                                         m_miJacIndex[ iTemp] = i;
02534                                         m_mdJacValue[ iTemp] = value[j];
02535                                         m_miJacStart[ index[j]]++;      
02536                                 }               
02537                         }                       
02538                 } 
02539         }
02540         //
02541         std::map<int, int>::iterator posVarIdx;
02542         // m_miJacStart[ i] is now equal to the correct m_miJacStart[ i] + m_miJacNumConTerms[ i], so readjust
02543         for (i = 0; i < iNumRowStarts - 1; i++ ){
02544                 m_miJacStart[ i] = m_miJacStart [ i] - m_miJacNumConTerms[ i] ; 
02545                 iTemp = m_miJacStart[ i] + m_miJacNumConTerms[ i];
02546                 // if the row is in the list of expression trees read in idicies and values
02547                 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
02548                         for(posVarIdx = (*m_mapExpressionTreesMod[ i]->mapVarIdx).begin(); posVarIdx 
02549                         != (*m_mapExpressionTreesMod[ i]->mapVarIdx).end(); ++posVarIdx){
02550                                 m_miJacIndex[ iTemp] = posVarIdx->first;
02551                                 m_mdJacValue[ iTemp] = 0;
02552                                 iTemp++;
02553                         }
02554                 }
02555         }
02556         #ifdef DEBUG
02557         std::cout << "HERE ARE ROW STARTS:" << std::endl;
02558         for (i = 0; i < iNumRowStarts; i++ ){
02559                 std::cout <<  m_miJacStart[ i] << "  "; 
02560         }
02561         std::cout << std::endl << std::endl;
02562         std::cout << "HERE ARE VARIABLE INDICIES:" << std::endl;
02563         for (i = 0; i < m_miJacStart[ iNumRowStarts - 1]; i++ ){
02564                 std::cout <<  m_miJacIndex[ i] << "  "; 
02565         }
02566         std::cout << std::endl << std::endl;
02567         std::cout << "HERE ARE VALUES:" << std::endl;
02568         for (i = 0; i < m_miJacStart[ iNumRowStarts - 1]; i++ ){
02569                 std::cout <<  m_mdJacValue[ i] << "  "; 
02570         }
02571         std::cout << std::endl << std::endl;
02572 
02573         std::cout << "HERE ARE NUMBER OF CONSTANT TERMS:" << std::endl;
02574         for (i = 0; i < iNumRowStarts - 1; i++ ){
02575                 std::cout <<  m_miJacNumConTerms[ i ] << "  ";  
02576         }
02577         std::cout << std::endl << std::endl;
02578         #endif
02579         return true;
02580 }//getSparseJacobianFromColumnMajor
02581 
02582 
02583 
02584 bool OSInstance::getSparseJacobianFromRowMajor( ){
02585         // we assume row major matrix
02586         if( m_bColumnMajor == true) return false;
02587         int iNumJacRowStarts = getConstraintNumber() + 1;
02588         std::map<int, int>::iterator posVarIdx; 
02589         int i,j, k;
02590         int *start = NULL;
02591         int *index = NULL;
02592         double *value = NULL;
02593         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02594                 start = this->instanceData->linearConstraintCoefficients->start->el;
02595                 index = this->instanceData->linearConstraintCoefficients->colIdx->el;
02596                 value = this->instanceData->linearConstraintCoefficients->value->el;
02597         }
02598         m_miJacStart = new int[ iNumJacRowStarts];
02599         m_miJacNumConTerms = new int[ getConstraintNumber()];
02600         OSnLNodePlus *nlNodePlus;
02601         OSnLNodeVariable *nlNodeVariable;
02602         //OSExpressionTree *expTree = NULL;
02603         // now initialize starts and variable index maps 
02604         for ( i = 0; i < iNumJacRowStarts; i++){                        
02605                 m_miJacStart [ i ] = 0;
02606                 // map the variables  in the nonlinear rows
02607                 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
02608                         // the following is equivalent to  m_treeRoot->getVariableIndexMap( i);
02609                         m_mapExpressionTreesMod[ i]->getVariableIndiciesMap();
02610                         
02611                 }
02612         }
02613         int loopLimit =  getConstraintNumber();
02614         // only execute the following code if there are linear constraint coefficients
02615         // determine the number of terms in constraint with constant partial derivative
02616         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02617                 // i is indexing rows (constrains) and j is indexing column numbers 
02618                 for (i = 0; i < loopLimit; i++){
02619                         m_miJacNumConTerms[ i] = 0;
02620                         for (j = start[i]; j < start[ i + 1 ]; j++){
02621                                 // determine if variable index[j] appears in the Expression Tree for row i
02622                                 // if we pass if test below then variable i is in the expresssion tree and we add
02623                                 // the linear term to the expession tree
02624                                 if( (m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ) &&
02625                                         ( (*m_mapExpressionTreesMod[ i]->mapVarIdx).find( index[ j]) != (*m_mapExpressionTreesMod[ i]->mapVarIdx).end()) ){
02626                                         // variable index[ j] appears in the expression tree for row i
02627                                         // add the coefficient corresponding to variable index[j] in row i to the expression tree       
02628                                         // define a new OSnLVariable and OSnLnodePlus 
02629                                         nlNodeVariable = new OSnLNodeVariable();
02630                                         nlNodeVariable->coef = value[ j];
02631                                         nlNodeVariable->idx = index[ j];
02632                                         nlNodePlus = new OSnLNodePlus();
02633                                         nlNodePlus->m_mChildren[ 0] = m_mapExpressionTreesMod[ i ]->m_treeRoot;
02634                                         nlNodePlus->m_mChildren[ 1] = nlNodeVariable;
02635                                         //expTree = new OSExpressionTree();
02636                                         //expTree->m_treeRoot = nlNodePlus ;
02637                                         //expTree->mapVarIdx = m_mapExpressionTreesMod[ i]->mapVarIdx;
02638                                         //m_mapExpressionTreesMod[ i ]  = expTree;      
02639                                         m_mapExpressionTreesMod[ i ]->m_treeRoot = nlNodePlus;
02640                                 }
02641                                 else{ 
02642                                         //the partial derivative of variable j in row i is constant
02643                                         m_miJacNumConTerms[ i]++;
02644                                 }                               
02645                         }
02646                 }
02647         }
02648         //
02649         m_miJacStart[0] = 0;
02650         for (i = 1; i < iNumJacRowStarts; i++ ){
02651                 if( m_mapExpressionTreesMod.find( i - 1) != m_mapExpressionTreesMod.end() ){
02652                         m_miJacStart[i] = m_miJacStart[i - 1] + (m_miJacNumConTerms[ i - 1] + (*m_mapExpressionTreesMod[ i - 1]->mapVarIdx).size() );
02653                 }
02654                 else{
02655                         m_miJacStart[i] = m_miJacStart[i - 1] + m_miJacNumConTerms[ i - 1];
02656                 }       
02657         }
02658         // we know how many constant terms and size of arrays
02659         // dimension miIndex and mdValue here
02660         m_iJacValueSize =       m_miJacStart[ iNumJacRowStarts - 1];
02661         m_miJacIndex = new int[  m_iJacValueSize];
02662         m_mdJacValue = new double[ m_iJacValueSize ];
02663         // now loop again and put in values and indicies
02664         // first put in the constant terms
02665         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02666                 for (i = 0; i < loopLimit; i++){
02667                         k = 0;
02668                         for (j = start[i]; j < start[ i + 1 ]; j++){
02669                                 if( (m_mapExpressionTreesMod.find( i) == m_mapExpressionTreesMod.end() ) ||
02670                                         ( (*m_mapExpressionTreesMod[ i]->mapVarIdx).find( index[ j]) == (*m_mapExpressionTreesMod[ i]->mapVarIdx).end()) ){
02671                                                 m_miJacIndex[ m_miJacStart[i] + k ] = index[ j];
02672                                                 m_mdJacValue[ m_miJacStart[i] + k ] = value[ j];
02673                                                 k++;
02674                                         }                       
02675                         }
02676                 }
02677         }
02678         // put in terms from the modified nonlinear expression tree
02679         for (i = 0; i < loopLimit; i++ ){
02680                 k = m_miJacStart[i] + m_miJacNumConTerms[ i ];
02681                 // if the row is in the list of expression trees read in idicies and values
02682                 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
02683                         for(posVarIdx = (*m_mapExpressionTreesMod[ i]->mapVarIdx).begin(); posVarIdx 
02684                         != (*m_mapExpressionTreesMod[ i]->mapVarIdx).end(); ++posVarIdx){
02685                                 m_miJacIndex[ k] = posVarIdx->first;
02686                                 m_mdJacValue[ k] = 0;
02687                                 k++;
02688                         }
02689                 }
02690         }
02691         #ifdef DEBUG
02692         std::cout << "HERE ARE ROW STARTS:" << std::endl;
02693         for (i = 0; i < iNumJacRowStarts; i++ ){
02694                 std::cout <<  m_miJacStart[ i] << "  "; 
02695         }
02696         std::cout << std::endl << std::endl;
02697         std::cout << "HERE ARE VARIABLE INDICIES:" << std::endl;
02698         for (i = 0; i < m_miJacStart[ iNumJacRowStarts - 1]; i++ ){
02699                 std::cout <<  m_miJacIndex[ i] << "  "; 
02700         }
02701         std::cout << std::endl << std::endl;
02702         std::cout << "HERE ARE VALUES:" << std::endl;
02703         for (i = 0; i < m_miJacStart[ iNumJacRowStarts - 1]; i++ ){
02704                 std::cout <<  m_mdJacValue[ i] << "  "; 
02705         }
02706         std::cout << std::endl << std::endl;
02707 
02708         std::cout << "HERE ARE NUMBER OF CONSTANT TERMS:" << std::endl;
02709         for (i = 0; i < iNumJacRowStarts - 1; i++ ){
02710                 std::cout <<  m_miJacNumConTerms[ i ] << "  ";  
02711         }
02712         std::cout << std::endl << std::endl;
02713         #endif
02714         return true;
02715 }//getSparseJacobianFromRowMajor
02716 
02717 OSExpressionTree* OSInstance::getLagrangianExpTree( ){
02718         if( m_bLagrangianExpTreeCreated == true) return m_LagrangianExpTree;
02719         // we calculate the Lagrangian for all the objectives and constraints
02720         // with nonlinear terms
02721         // first initialize everything for nonlinear work
02722         if(m_bSparseJacobianCalculated == false) getJacobianSparsityPattern( ); 
02723         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02724         OSnLNodeTimes* nlNodeTimes = NULL;
02725         OSnLNodeVariable* nlNodeVariable = NULL;
02726         OSnLNodeSum* nlNodeSum = NULL;
02727         int numChildren = 0;
02728         int rowIdx;
02729         // create the sum node
02730         nlNodeSum = new OSnLNodeSum();
02731         nlNodeSum->inumberOfChildren = m_mapExpressionTreesMod.size();
02732         nlNodeSum->m_mChildren = new OSnLNode*[ nlNodeSum->inumberOfChildren];
02733         // create and expression tree for the sum node
02734         m_LagrangianExpTree = new OSExpressionTree();
02735         m_LagrangianExpTree->m_treeRoot = nlNodeSum;
02736         // now create the children of the sum node
02737         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02738                 // this variable is the Lagrange multiplier
02739                 nlNodeVariable = new OSnLNodeVariable();
02740                 nlNodeVariable->coef = 1.;
02741                 // get the correct index --
02742                 // for rowIdx = 0, ..., m - 1 set idx = numVar + rowIdx
02743                 rowIdx = posMapExpTree->first;
02744                 if(rowIdx >= 0){
02745                         nlNodeVariable->idx = instanceData->variables->numberOfVariables + rowIdx;
02746                 }
02747                 else{
02748                         nlNodeVariable->idx = instanceData->variables->numberOfVariables + 
02749                         instanceData->constraints->numberOfConstraints + (abs(rowIdx) - 1);
02750                 }
02751                 // now create a times multiply the new variable times the root of the expression tree
02752                 nlNodeTimes = new OSnLNodeTimes();
02753                 nlNodeTimes->m_mChildren[ 0] = nlNodeVariable;
02754                 nlNodeTimes->m_mChildren[ 1] = m_mapExpressionTreesMod[ posMapExpTree->first ]->m_treeRoot;     
02755                 // the times node is the new child
02756                 nlNodeSum->m_mChildren[ numChildren] = nlNodeTimes;
02757                 numChildren++;
02758         }       
02759         // get a variable index map for the expression tree
02760         m_LagrangianExpTree->getVariableIndiciesMap();
02761         // print out the XML for this puppy
02762         //std::cout << "Here comes the Lagrangian Tree" << std::endl;
02763         //std::cout << m_LagrangianExpTree->m_treeRoot->getNonlinearExpressionInXML() << std::endl;
02764         //
02765         m_bLagrangianExpTreeCreated = true;
02766         return m_LagrangianExpTree;
02767 }//getLagrangianExpTree
02768 
02769 std::map<int, int> OSInstance::getAllNonlinearVariablesIndexMap( ){
02770         if(m_bAllNonlinearVariablesIndex == true) return m_mapAllNonlinearVariablesIndex;
02771         //loop over the map of expression tree and get a unique listing of all variables
02772         // put these in the map m_mapAllNonlinearVariablesIndex
02773         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02774         std::map<int, int>::iterator posVarIdx;
02775         OSExpressionTree *expTree;
02776         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02777                 // get the index map for the expression tree
02778                 expTree = posMapExpTree->second;
02779                 if(expTree->m_bIndexMapGenerated == false)expTree->getVariableIndiciesMap();
02780                 for(posVarIdx = (*expTree->mapVarIdx).begin(); posVarIdx != (*expTree->mapVarIdx).end(); ++posVarIdx){
02781                         if( m_mapAllNonlinearVariablesIndex.find( posVarIdx->first) == m_mapAllNonlinearVariablesIndex.end() ){
02782                         // add the variable to the Lagragian map
02783                         m_mapAllNonlinearVariablesIndex[ posVarIdx->first] = 1;
02784                         }
02785                 }
02786         }
02787         m_miNonLinearVarsReverseMap = new int[m_mapAllNonlinearVariablesIndex.size()];
02788         // now order appropriately
02789         int kount = 0;
02790         //std::cout << "HERE IS THE LAGRANGIANN VARIABLE MAPPING" << std::endl;
02791         for(posVarIdx = m_mapAllNonlinearVariablesIndex.begin(); posVarIdx !=m_mapAllNonlinearVariablesIndex.end(); ++posVarIdx){
02792                 posVarIdx->second = kount;
02793                 m_miNonLinearVarsReverseMap[ kount] = posVarIdx->first;
02794                 kount++;
02795                 #ifdef DEBUG
02796                 std::cout <<  "POSITION FIRST =  "  << posVarIdx->first ;
02797                 std::cout <<  "   POSITION SECOND = "  << posVarIdx->second << std::endl;
02798                 #endif
02799         }
02800         m_iNumberOfNonlinearVariables = kount;
02801         //std::cout <<  "NUMBER OF NONLINEAR VARIABLES =  "  << kount ;
02802         m_bAllNonlinearVariablesIndex = true;
02803         return m_mapAllNonlinearVariablesIndex;
02804 }//getAllNonlinearVariablesIndexMap     
02805 
02806 SparseHessianMatrix* OSInstance::getLagrangianHessianSparsityPattern( ){
02807         // fill in the nonzeros in the sparse Hessian
02808         if( m_bLagrangianSparseHessianCreated == true) return m_LagrangianSparseHessian;
02809         if( m_iNumberOfNonlinearVariables == 0) return NULL;
02810         if( m_binitForAlgDiff == false ) initForAlgDiff();
02811         unsigned int i = 0;
02812         int numNonz = 0;
02813         // Create the CppAD function if necessary
02814         //
02815         std::vector<double> vx;
02816         std::map<int, int>::iterator posMap1, posMap2;  
02817         if( (m_bCppADFunIsCreated == false || m_bCppADMustReTape == true )  && (m_mapExpressionTreesMod.size() > 0) ) {
02818                 for(posMap1 = m_mapAllNonlinearVariablesIndex.begin(); posMap1 != m_mapAllNonlinearVariablesIndex.end(); ++posMap1){
02819                         vx.push_back( 1.0) ;
02820                 }
02821                 createCppADFun( vx);
02822         }
02823         //
02824         // Use CppAD to do a forward sparsity calculation
02825         std::vector<bool> r(m_iNumberOfNonlinearVariables * m_iNumberOfNonlinearVariables);
02826         unsigned int j;
02827         for(i = 0; i < m_iNumberOfNonlinearVariables; i++) { 
02828                 for(j = 0; j < m_iNumberOfNonlinearVariables; j++)
02829                         r[ i * m_iNumberOfNonlinearVariables + j ] = false;
02830                         r[ i * m_iNumberOfNonlinearVariables + i] = true;
02831         }
02832         // compute sparsity pattern for J(x) = F^{(1)} (x)
02833         (*Fad).ForSparseJac(m_iNumberOfNonlinearVariables, r);
02834         //
02835         //now the second derivative
02836         unsigned int m = m_mapExpressionTreesMod.size();
02837         std::vector<bool> e( m);
02838         //Vector s(m);
02839         for(i = 0; i < m; i++) e[i] = true;
02840         std::cout << "Computing Sparse Hessian" << std::endl;
02841         //m_vbLagHessNonz holds the sparsity pattern Lagrangian of the Hessian
02842         m_vbLagHessNonz = (*Fad).RevSparseHes(m_iNumberOfNonlinearVariables, e);
02843         for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
02844                 //std::cout << "Row " << i << "  of Hessian " << std::endl;
02845                 for(j = i; j < m_iNumberOfNonlinearVariables; j++){
02846                         if(m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j]  == true) numNonz++;
02847                         //std::cout << m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j] <<  "  " ;
02848                 }
02849                 //std::cout << std::endl;
02850         }
02851         //std::cout << "Lagrangian Hessian Nonzeros = " << numNonz << std::endl;
02852         i = 0;
02853         // now that we have the dimension create SparseHessianMatrix (upper triangular)
02854         m_LagrangianSparseHessian = new SparseHessianMatrix();
02855         m_LagrangianSparseHessian->bDeleteArrays = true;
02856         m_LagrangianSparseHessian->hessDimension = numNonz;
02857         m_LagrangianSparseHessian->hessRowIdx = new int[m_LagrangianSparseHessian->hessDimension];
02858         m_LagrangianSparseHessian->hessColIdx = new int[m_LagrangianSparseHessian->hessDimension];
02859         m_LagrangianSparseHessian->hessValues = new double[m_LagrangianSparseHessian->hessDimension];
02860         //std::cout << "HESSIAN DIMENSION = " << m_LagrangianSparseHessian->hessDimension << std::endl;
02861         numNonz = 0;
02862         for(posMap1 = m_mapAllNonlinearVariablesIndex.begin(); posMap1 != m_mapAllNonlinearVariablesIndex.end(); ++posMap1){
02863                 //std::cout << "posMap1->first  " << posMap1->first << std::endl;
02864                 j = i;
02865                 for(posMap2 = posMap1; posMap2 != m_mapAllNonlinearVariablesIndex.end(); ++posMap2){
02866                         //std::cout << "posMap2->first  " << posMap2->first << std::endl;
02867                         if(m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j] == true){
02868                                 *(m_LagrangianSparseHessian->hessRowIdx + numNonz) = posMap1->first;
02869                                 *(m_LagrangianSparseHessian->hessColIdx + numNonz) = posMap2->first;
02870                                 numNonz++;
02871                         }
02872                         //std::cout << m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j] <<  "  " << std::endl;
02873                         j++;
02874                 }
02875                 i++;
02876         }
02877         #ifdef DEBUG
02878         std::cout << "HESSIAN SPARSITY PATTERN" << std::endl;
02879         int kj;
02880         for(kj = 0; kj < m_LagrangianSparseHessian->hessDimension; kj++){
02881                 std::cout <<  "Row Index = " << *(m_LagrangianSparseHessian->hessRowIdx + kj) << std::endl;
02882                 std::cout <<  "Column Index = " << *(m_LagrangianSparseHessian->hessColIdx + kj) << std::endl;
02883         }
02884         #endif
02885         //
02886         m_bLagrangianSparseHessianCreated = true;
02887         return m_LagrangianSparseHessian;
02888 }//getLagrangianHessianSparsityPattern
02889 
02890 
02891 void OSInstance::duplicateExpressionTreesMap(){
02892         //std::cout << "I AM IN DUPLICATE EXPRSSION TREES MAP" << std::endl;
02893         // we do this so that we can keep the integrity of the original formulation
02894         if(m_bDuplicateExpressionTreesMap == false){ 
02895                 // first make sure the map was created
02896                 if( m_bProcessExpressionTrees == false) getAllNonlinearExpressionTrees();
02897                 m_mapExpressionTreesMod = m_mapExpressionTrees;
02898                 m_bDuplicateExpressionTreesMap = true;
02899                 return;
02900         }
02901         else{
02902                 return;
02903         }
02904 }//duplicateExpressionTreesMap
02905 
02906 
02907 bool OSInstance::createCppADFun(std::vector<double> vdX){
02908         if(m_bCppADFunIsCreated == true) return true;
02909         //if( m_bNonLinearStructuresInitialized == false) initializeNonLinearStructures( );
02910         if(m_binitForAlgDiff == false) initForAlgDiff();
02911         
02912         //if( m_bAllNonlinearVariablesIndex == false) getAllNonlinearVariablesIndexMap( );
02913         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02914         unsigned int i;
02915         size_t n = vdX.size();
02916         // declare a CppAD vector and fill it in
02917         CppADvector< AD<double> > vdaX( n );
02918         for(i = 0; i < n; i++){
02919                 vdaX[ i] = vdX[ i];
02920                 //std::cout << "vdX =  " << vdX[ i] << std::endl;
02921         }
02922         // declare the independent variables and start recording
02923         CppAD::Independent( vdaX);
02929         CppAD::vector< AD<double> > m_vFG;      
02930         int kount = 0;
02931         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){  
02932                 m_vFG.push_back( (posMapExpTree->second)->m_treeRoot->constructCppADTape(&m_mapAllNonlinearVariablesIndex, &vdaX) );
02933                 //std::cout << "PUSHING BACK EXPRESSION NUMBER " << posMapExpTree->first << std::endl;
02934                 if( m_mapCppADFunRangeIndex.find( posMapExpTree->first) == m_mapCppADFunRangeIndex.end() ){
02935                         // count which nonlinear obj/constraint this is
02936                         m_mapCppADFunRangeIndex[ posMapExpTree->first] = kount;
02937                         kount++;
02938                 }
02939         }       
02940         //create the function and stop recording
02941         std::cout << "create the function and stop recording"  << std::endl;
02942         Fad = new CppAD::ADFun<double>(vdaX, m_vFG);
02943         std::cout << "range space dimension =  " << m_vFG.size() << std::endl;
02944         // no forward sweeps done yet
02945         m_iHighestTaylorCoeffOrder = -1;
02946         m_bCppADFunIsCreated = true;
02947         return true;
02948 }//end createCppADFun
02949 
02950 
02951 std::vector<double> OSInstance::forwardAD(int p, std::vector<double> vdX){
02952         try{
02953                 // make sure a CppADFun has been created
02954                 if(m_bCppADFunIsCreated == false) createCppADFun( vdX);
02955                 if(p > (m_iHighestTaylorCoeffOrder + 1) ) throw 
02956                         ErrorClass( "trying to calculate a p order forward when p-1 Taylor coefficient not available");
02957                 // adjust the order of the Taylor coefficient
02958                 m_iHighestTaylorCoeffOrder = p; 
02959                 m_iHighestOrderEvaluated = p;
02960                 //for(int i  = 0; i < vdX.size(); i++){
02961                         //std::cout << "ForwardAD Primal Variables " << i   << " " << vdX[ i] << std::endl;
02962                 //}
02963                 return (*Fad).Forward(p, vdX);
02964         }
02965         catch(const ErrorClass& eclass){
02966                 throw ErrorClass( eclass.errormsg);
02967         }  
02968 }//end forwardAD
02969 
02970 
02971 std::vector<double> OSInstance::reverseAD(int p, std::vector<double> vdlambda){
02972         try{
02973                 if(p == 0) throw 
02974                         ErrorClass( "reverseAD must have p >= 1");
02975                 if(p > (m_iHighestTaylorCoeffOrder + 1) ) throw 
02976                         ErrorClass( "trying to calculate a p order reverse when p-1 Taylor coefficient not available");
02977                 //for(int i  = 0; i < vdlambda.size(); i++){
02978                 //      std::cout << "ReverseAD Multiplier " << i   << " " << vdlambda[ i] << std::endl;
02979                 //}
02980                 m_iHighestOrderEvaluated = p;
02981                 return (*Fad).Reverse(p, vdlambda);
02982         }
02983         catch(const ErrorClass& eclass){
02984                 throw ErrorClass( eclass.errormsg);
02985         }  
02986 }//end forwardAD
02987 
02988 bool OSInstance::getIterateResults( double *x, double *objLambda, double* conMultipliers, 
02989                 bool new_x, int highestOrder){
02990         try{ 
02991                 if( m_binitForAlgDiff == false) initForAlgDiff();
02992                 std::map<int, int>::iterator posVarIndexMap;
02993                 
02994                 if(new_x == true){
02995                         if( m_vdX.size() > 0) m_vdX.clear();
02996                         for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
02997                                 m_vdX.push_back( x[ posVarIndexMap->first]) ;
02998                         }
02999                         if( (m_bCppADFunIsCreated == false || m_bCppADMustReTape == true )  && (m_mapExpressionTreesMod.size() > 0) ) {
03000                                 createCppADFun( m_vdX);
03001                         }       
03002                 }
03003                 else{ // make sure vector not empty -- this could happen if we have linear obj and nonlinear constraints
03004                         if( m_vdX.size() == 0) {
03005                                 for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
03006                                         m_vdX.push_back( x[ posVarIndexMap->first]) ;
03007                                 }
03008                                 if( (m_bCppADFunIsCreated == false || m_bCppADMustReTape == true )  && (m_mapExpressionTreesMod.size() > 0) ) {
03009                                         createCppADFun( m_vdX);
03010                                 }       
03011                         }
03012                 }       
03013                 switch( highestOrder){          
03014                         case 0: 
03015                                 if(new_x == true || m_iHighestOrderEvaluated < 0){      
03016                                         if(bUseExpTreeForFunEval == true){
03017                                                 calculateAllConstraintFunctionValues( x, new_x);
03018                                                 calculateAllObjectiveFunctionValues( x, new_x);
03019                                         }
03020                                         else{
03021                                                 getZeroOrderResults(x, objLambda, conMultipliers);
03022                                         }
03023 
03024                                 }
03025                                 break;  
03026                         case 1:
03027                                 if(new_x == true || m_iHighestOrderEvaluated < 0)       
03028                                         getZeroOrderResults(x, objLambda, conMultipliers);
03029                                 if(new_x == true || m_iHighestOrderEvaluated < 1)       
03030                                         getFirstOrderResults(x, objLambda, conMultipliers);
03031                                 break;
03032                         case 2: 
03033                                 if(new_x == true || m_iHighestOrderEvaluated < 0)       
03034                                         getZeroOrderResults(x, objLambda, conMultipliers);
03035                                 if(new_x == true || m_iHighestOrderEvaluated < 2)       
03036                                         getSecondOrderResults(x, objLambda, conMultipliers);
03037                                 break;
03038                         default:
03039                                 throw ErrorClass("Derivative should be order 0, 1, or 2");      
03040                 }//end switch
03041                 return true;
03042         }
03043         catch(const ErrorClass& eclass){
03044                 throw ErrorClass( eclass.errormsg);
03045         }  
03046 }//end getIterateResults
03047 
03048 
03049 bool OSInstance::getZeroOrderResults(double *x, double *objLambda, double *conMultipliers){
03050         try{ 
03051                 // initialize everything
03052                 int i, j, rowNum, objNum;
03053                 if( m_mapExpressionTreesMod.size() > 0){
03054                         m_vdYval = this->forwardAD(0, m_vdX);   
03055                 }
03056                 // now get all function and constraint values using forward result
03057                 for(rowNum = 0; rowNum < m_iConstraintNumber; rowNum++){
03058                         m_mdConstraintFunctionValues[ rowNum] = 0.0;
03059                         if( m_mapExpressionTreesMod.find( rowNum) != m_mapExpressionTreesMod.end() ){
03060                                 m_mdConstraintFunctionValues[ rowNum] = m_vdYval[  m_mapCppADFunRangeIndex[ rowNum]];
03061                         }
03062                         // now the linear part
03063                         // be careful, loop over only the constant terms in sparseJacMatrix
03064                         i = m_sparseJacMatrix->starts[ rowNum];
03065                         j = m_sparseJacMatrix->starts[ rowNum + 1 ];
03066                         while ( (i - m_sparseJacMatrix->starts[ rowNum])  < m_sparseJacMatrix->conVals[ rowNum] ){
03067                                 m_mdConstraintFunctionValues[ rowNum] += m_sparseJacMatrix->values[ i]*x[ m_sparseJacMatrix->indexes[ i] ];
03068                                 i++;
03069                         }       
03070                         // add in the constraint function constant
03071                         m_mdConstraintFunctionValues[ rowNum] += m_mdConstraintConstants[ rowNum ];
03072                         #ifdef DEBUG
03073                         std::cout << "Constraint " <<  rowNum << " function value =  " << m_mdConstraintFunctionValues[ rowNum ] << std::endl;
03074                         #endif
03075                 }
03076                 // now get the objective function values from the forward result
03077                 for(objNum = 0; objNum < m_iObjectiveNumber; objNum++){
03078                         m_mdObjectiveFunctionValues[ objNum] = 0.0;
03079                         if( m_mapExpressionTreesMod.find( -objNum -1) != m_mapExpressionTreesMod.end() ){
03080                                 m_mdObjectiveFunctionValues[ objNum] = m_vdYval[ objNum];
03081                         }
03082                         for(i = 0; i < m_iVariableNumber; i++){
03083                                 m_mdObjectiveFunctionValues[ objNum] += m_mmdDenseObjectiveCoefficients[ objNum][i]*x[ i];
03084                         }
03085                         #ifdef DEBUG
03086                         std::cout << "Objective " << objNum << " function value =  " << m_mdObjectiveFunctionValues[ objNum] << std::endl;
03087                         #endif
03088                 }
03089         return true;
03090         }//end try
03091         catch(const ErrorClass& eclass){
03092                 throw ErrorClass( eclass.errormsg);
03093         }  
03094 }//end getZeroOrderResults
03095 
03096 
03097 
03098 bool OSInstance::getFirstOrderResults(double *x, double *objLambda, double *conMultipliers){
03099         try{
03100                 // initialize everything
03101                 unsigned int i, j;
03102                 int rowNum,  jacIndex;
03103                 unsigned int jstart, jend;
03104                 int idx;
03105                 OSExpressionTree *expTree = NULL;
03106                 int domainIdx = 0;      
03107                 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03108                 std::map<int, int>::iterator posVarIdx;
03109                         
03115                 if(m_iNumberOfNonlinearVariables >= m_mapExpressionTreesMod.size() ){
03116                         // calculate the gradient by doing a reverse sweep over each row
03117                         // loop over the constraints that have a nonlinear term and get their gradients
03118                         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03119                                 idx = posMapExpTree->first;
03120                                 // we are considering only constraints, not objective function
03121                                 if(idx >= 0){
03122                                         m_vdRangeUnitVec[ domainIdx] = 1.;
03123                                         m_mapExpressionTreesMod[ idx]->getVariableIndiciesMap(); 
03124                                         m_vdYjacval = this->reverseAD(1, m_vdRangeUnitVec);
03125                                         // check size
03126                                         jstart = m_miJacStart[ idx] + m_miJacNumConTerms[ idx];
03127                                         jend = m_miJacStart[ idx + 1 ];
03128                                         if( (*m_mapExpressionTreesMod[ idx]->mapVarIdx).size() != (jend - jstart)) throw 
03129                                         ErrorClass("number of partials not consistent");
03130                                         j = 0;
03131                                         jacIndex = 0;
03132                                         for(posVarIdx = m_mapAllNonlinearVariablesIndex.begin(); posVarIdx 
03133                                                 != m_mapAllNonlinearVariablesIndex.end(); ++posVarIdx){
03134                                                 //std::cout << "Constraint Function Jacobian Values" << "For Constraint  " << idx  << std::endl;
03135                                                 //std::cout << "Jac Val for index " << posVarIdx->first  << " = " << m_vdYjacval[ jacIndex] << std::endl;
03136                                                 //if(m_miJacIndex[ jstart] != posVarIdx->first) throw ErrorClass("error calculating Jacobian matrix");
03137                                                 // we are working with variable posVarIdx->first in the original variable space
03138                                                 // we need to see which variable this is in the individual constraint map
03139                                                 if( (*m_mapExpressionTreesMod[ idx]->mapVarIdx).find( posVarIdx->first) != (*m_mapExpressionTreesMod[ idx]->mapVarIdx).end()){
03140                                                         m_mdJacValue[ jstart] = m_vdYjacval[ jacIndex];
03141                                                         jstart++;
03142                                                         j++;
03143                                                 }
03144                                                 jacIndex++;
03145                                         }
03146                                         
03147                                         m_vdRangeUnitVec[ domainIdx] = 0.;
03148                                         domainIdx++;
03149                                 }
03150                                 else{    // we have an objective function
03151 
03152                                         domainIdx++;
03153                                 }
03154                         }
03155                 }
03156                 else{  
03157                         // calculate the gradients using a forward sweep over all the variables.                
03158                         for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
03159                                 m_vdDomainUnitVec[i] = 1.;     
03160                                 rowNum = 0;
03161                                 if( m_mapExpressionTreesMod.size() > 0){          
03162                                         m_vdYjacval = this->forwardAD(1, m_vdDomainUnitVec);
03163                                 } 
03164                                 // fill in Jacobian here, we have column i 
03165                                 // start Jacobian calculation
03166                                 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03167                                         idx = posMapExpTree->first;
03168                                         // we are considering only constraints, not objective function
03169                                         if(idx >= 0){
03170                                                 //figure out original variable this corresponds to
03171                                                 //then use (*m_mapExpressionTreesMod[ idx]->mapVarIdx) to figure out which variable it is within row idx
03172                                                 //m_mapAllNonlinearVariablesIndex
03173                                                 //std::cout << "This is the following variable in the expression tree  " <<  (*m_mapExpressionTreesMod[ idx]->mapVarIdx)[ m_miNonLinearVarsReverseMap[ i]]<< std::endl;                         
03174                                                 expTree = m_mapExpressionTreesMod[ idx];                
03175                                                 if( (*expTree->mapVarIdx).find( m_miNonLinearVarsReverseMap[ i]) != (*expTree->mapVarIdx).end()  ){             
03176                                                         jacIndex = (*m_mapExpressionTreesMod[ idx]->mapVarIdx)[ m_miNonLinearVarsReverseMap[ i]];
03177                                                         jstart = m_miJacStart[ idx] + m_miJacNumConTerms[ idx];
03178                                                         // kipp change 1 to number of objective functions
03179                                                         m_mdJacValue[ jstart + jacIndex] = m_vdYjacval[m_iObjectiveNumberNonlinear + rowNum];
03180                                                 }
03181                                                 rowNum++;
03182                                         }//end Jacobian calculation
03183 
03184                         }                       
03185                         //
03186                         m_vdDomainUnitVec[i] = 0.;
03187                         }
03188                 }
03189                 #ifdef DEBUG
03190                 int k;
03191                 std::cout  << "JACOBIAN DATA " << std::endl;
03192                 for(idx = 0; idx < m_iConstraintNumber; idx++){
03193                         for(k = *(m_sparseJacMatrix->starts + idx); k < *(m_sparseJacMatrix->starts + idx + 1); k++){
03194                                 std::cout << "row idx = " << idx <<  "  col idx = "<< *(m_sparseJacMatrix->indexes + k)
03195                                 << " value = " << *(m_sparseJacMatrix->values + k) << std::endl;
03196                         }
03197                 }
03198 
03199                 #endif
03200                 return true;
03201         }//end try
03202         catch(const ErrorClass& eclass){
03203                 throw ErrorClass( eclass.errormsg);
03204         } 
03205 }// end getFirstOrderResults
03206                         
03207 
03208 bool OSInstance::getSecondOrderResults(double *x, double *objLambda, double *conMultipliers){
03209         try{
03210                 // initialize everything
03211                 unsigned int i, j;
03212                 int rowNum,  jacIndex;
03213                 int jstart,  idx;
03214                 OSExpressionTree *expTree = NULL;
03215                 int hessValuesIdx = 0;  
03216                 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03217                 std::map<int, int>::iterator posVarIndexMap;
03218                 if( conMultipliers == NULL) throw ErrorClass("cannot have a null vector of lagrange multipliers when calculating Hessian of Lagrangian");
03219                 if( m_vdLambda.size() > 0) m_vdLambda.clear();
03220                 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){  
03221                         if( posMapExpTree->first >= 0){
03222                                 m_vdLambda.push_back( conMultipliers[ posMapExpTree->first]);
03223                         }
03224                         else{
03225                                 // kipp correct when there is more than one obj
03226                                 m_vdLambda.push_back( objLambda[ abs(posMapExpTree->first) - 1] );
03227                         }
03228                 }
03229                 for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
03230                         m_vdDomainUnitVec[i] = 1.;     
03231                         rowNum = 0;
03232                         if( m_mapExpressionTreesMod.size() > 0){          
03233                                 m_vdYjacval = this->forwardAD(1, m_vdDomainUnitVec);
03234                         } 
03235                         // fill in Jacobian here, we have column i 
03236                         // start Jacobian calculation
03237                         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03238                                 idx = posMapExpTree->first;
03239                                 // we are considering only constraints, not objective function
03240                                 if(idx >= 0){
03241                                         //figure out original variable this corresponds to
03242                                         //then use (*m_mapExpressionTreesMod[ idx]->mapVarIdx) to figure out which variable it is within row idx
03243                                         //m_mapAllNonlinearVariablesIndex
03244                                         //std::cout << "This is the following variable in the expression tree  " <<  (*m_mapExpressionTreesMod[ idx]->mapVarIdx)[ m_miNonLinearVarsReverseMap[ i]]<< std::endl;                         
03245                                         expTree = m_mapExpressionTreesMod[ idx];                
03246                                         if( (*expTree->mapVarIdx).find( m_miNonLinearVarsReverseMap[ i]) != (*expTree->mapVarIdx).end()  ){             
03247                                                 jacIndex = (*m_mapExpressionTreesMod[ idx]->mapVarIdx)[ m_miNonLinearVarsReverseMap[ i]];
03248                                                 jstart = m_miJacStart[ idx] + m_miJacNumConTerms[ idx];
03249                                                 m_mdJacValue[ jstart + jacIndex] = m_vdYjacval[m_iObjectiveNumberNonlinear + rowNum];
03250                                         }
03251                                         rowNum++;
03252                                 }//end Jacobian calculation
03253                                 else{
03254                                         // see if we have the objective function of interest
03255                                         //kipp fix if more than one obj
03256                                                 m_mmdObjGradient[  (abs( idx) - 1)][ m_miNonLinearVarsReverseMap[ i]] = m_vdYjacval[ (abs( idx) - 1)] + 
03257                                                 m_mmdDenseObjectiveCoefficients[  (abs( idx) - 1)][ m_miNonLinearVarsReverseMap[ i]];                                   
03258                         }//end Obj gradient calculation 
03259                 }                       
03260                 // now calculate the Hessian
03261                 if( m_mapExpressionTreesMod.size() > 0){   
03262                         m_vdw = reverseAD(2, m_vdLambda);   // derivtative of partial
03263                 }
03264                 for(j = i; j < m_iNumberOfNonlinearVariables; j++){
03265                         if( m_vbLagHessNonz[i*m_iNumberOfNonlinearVariables + j] == true){
03266                                 m_LagrangianSparseHessian->hessValues[ hessValuesIdx] =  m_vdw[  j*2 + 1];
03267                                 #ifdef DEBUG
03268                                 std::cout << "reverse 2 " << m_LagrangianSparseHessian->hessValues[ hessValuesIdx] << std::endl;
03269                                 #endif
03270                                 hessValuesIdx++;
03271                         }
03272                 }
03273                 //
03274                 //
03275                 m_vdDomainUnitVec[i] = 0.;
03276         }
03277 
03278         #ifdef DEBUG
03279         int k;
03280         std::cout  << "JACOBIAN DATA " << std::endl;
03281         for(idx = 0; idx < m_iConstraintNumber; idx++){
03282                 for(k = *(m_sparseJacMatrix->starts + idx); k < *(m_sparseJacMatrix->starts + idx + 1); k++){
03283                         std::cout << "row idx = " << idx <<  "  col idx = "<< *(m_sparseJacMatrix->indexes + k)
03284                         << " value = " << *(m_sparseJacMatrix->values + k) << std::endl;
03285                 }
03286         }
03287         #endif
03288         return true;
03289         }//end try
03290         catch(const ErrorClass& eclass){
03291                 throw ErrorClass( eclass.errormsg);
03292         } 
03293 }// end getSecondOrderResults
03294 
03295 bool OSInstance::initForAlgDiff(){
03296         if( m_binitForAlgDiff == true ) return true;
03297         initializeNonLinearStructures( );
03298         initObjGradients();
03299         getAllNonlinearVariablesIndexMap( );
03300         if(m_bSparseJacobianCalculated  == false) getJacobianSparsityPattern();
03301         //see if we need to retape 
03302         //loop over expression tree and see if one requires it
03303         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03304         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03305                 if(posMapExpTree->second->bCppADMustReTape == true) m_bCppADMustReTape = true;
03306         }                               
03307 
03308         #ifdef DEBUG
03309         std::cout << "RETAPE ==  " << m_bCppADMustReTape << std::endl;
03310         #endif
03311         unsigned int i;
03312         for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
03313                 m_vdDomainUnitVec.push_back( 0.0 );
03314         }
03315         for(i = 0; i < m_mapExpressionTreesMod.size(); i++){
03316                 m_vdRangeUnitVec.push_back( 0.0 );
03317         }
03318         m_binitForAlgDiff = true;
03319         return true;
03320 }//end initForAlgDiff
03321 
03322 bool OSInstance::initObjGradients(){
03323         int i, j;
03324         int m, n;
03325         m = getObjectiveNumber();
03326         n = getVariableNumber();
03327         getDenseObjectiveCoefficients();
03328         m_mmdObjGradient = new double*[m];
03329         for(i = 0; i < m; i++){
03330                 m_mmdObjGradient[i] = new double[n];
03331                 for(j = 0; j < n; j++){
03332                         m_mmdObjGradient[i][j] =  m_mmdDenseObjectiveCoefficients[ i][j];
03333                         #ifdef DEBUG
03334                         std::cout << "m_mmdObjGradient[i][j] = " << m_mmdObjGradient[i][j]  << std::endl;
03335                         #endif
03336                 }
03337         }
03338         return true;
03339 }//end initObjGradients

Generated on Tue Sep 30 03:01:24 2008 by  doxygen 1.4.7