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

Go to the documentation of this file.
00001 /* $Id: OSInstance.cpp 2831 2009-07-21 07:34:44Z kmartin $ */
00020 #include "OSInstance.h"
00021 #include "OSMathUtil.h"
00022 #include "OSErrorClass.h"
00023 #include "OSParameters.h"
00024 
00025 #include<stack>
00026 #include<iostream>  
00027 #include<sstream>
00028 
00029 //#define DEBUG
00030  
00031 using namespace std;
00032 using std::ostringstream; 
00033 
00034 
00035 OSInstance::OSInstance(): 
00036         bVariablesModified(false),
00037         bObjectivesModified(false),
00038         bConstraintsModified(false),
00039         bAMatrixModified(false),
00040         m_sInstanceName(""),
00041         m_sInstanceSource(""),  
00042         m_sInstanceDescription(""),
00043         m_bProcessVariables(false),
00044         m_iVariableNumber(-1),
00045         m_iNumberOfIntegerVariables( 0),
00046         m_iNumberOfBinaryVariables( 0),
00047         m_iNumberOfStringVariables( 0),
00048         m_iNumberOfQuadraticRowIndexes( 0),
00049         m_bQuadraticRowIndexesProcessed( false),
00050         m_miQuadRowIndexes( NULL),
00051         m_iNumberOfNonlinearExpressionTreeIndexes( 0),
00052         m_bNonlinearExpressionTreeIndexesProcessed( false),
00053         m_miNonlinearExpressionTreeIndexes( NULL),
00054         m_iNumberOfNonlinearExpressionTreeModIndexes( 0),
00055         m_bNonlinearExpressionTreeModIndexesProcessed( false),
00056         m_miNonlinearExpressionTreeModIndexes( NULL),           
00057         m_msVariableNames(NULL),
00058         //m_mdVariableInitialValues(NULL), -- deprecated
00059         //m_msVariableInitialStringValues(NULL), -- deprecated
00060         m_mcVariableTypes(NULL),
00061         m_mdVariableLowerBounds(NULL),
00062         m_mdVariableUpperBounds(NULL),
00063         m_bProcessObjectives(false),
00064         m_iObjectiveNumber(-1),
00065         m_iObjectiveNumberNonlinear( 0),
00066         m_msObjectiveNames(NULL),
00067         m_msMaxOrMins(NULL),
00068         m_miNumberOfObjCoef(NULL),
00069         m_mdObjectiveConstants(NULL),
00070         m_mdObjectiveWeights(NULL),
00071         m_mObjectiveCoefficients(NULL),
00072         m_bGetDenseObjectives(false),
00073         m_mmdDenseObjectiveCoefficients(NULL),
00074         m_bProcessConstraints(false),
00075         m_iConstraintNumber(-1),
00076         m_iConstraintNumberNonlinear( 0),       
00077         m_msConstraintNames(NULL),
00078         m_mdConstraintLowerBounds(NULL),
00079         m_mdConstraintUpperBounds(NULL),
00080         m_mdConstraintConstants( NULL),
00081         m_mcConstraintTypes(NULL),
00082         m_bProcessLinearConstraintCoefficients(false),  
00083         m_iLinearConstraintCoefficientNumber(-1),
00084         m_bColumnMajor(true),
00085         m_binitForAlgDiff( false),      
00086         m_linearConstraintCoefficientsInColumnMajor(NULL),
00087         m_linearConstraintCoefficientsInRowMajor(NULL), 
00088         m_bProcessQuadraticTerms(false),
00089         m_iQuadraticTermNumber(-1),
00090         m_mdConstraintFunctionValues( NULL),
00091         m_mdObjectiveFunctionValues( NULL),
00092         m_iJacValueSize( 0),
00093         m_miJacStart( NULL),
00094         m_miJacIndex( NULL),
00095         m_mdJacValue( NULL),
00096         m_miJacNumConTerms( NULL),
00097         m_sparseJacMatrix( NULL),       
00098         m_iHighestTaylorCoeffOrder(-1), 
00099         m_quadraticTerms( NULL),        
00100         m_bQTermsAdded( false), 
00101         m_iNumberOfNonlinearVariables( 0),
00102         m_bProcessNonlinearExpressions( false),
00103         m_iNonlinearExpressionNumber( -1),              
00104         m_miNonlinearExpressionIndexes( NULL),
00105         m_bProcessExpressionTrees( false),
00106         m_bProcessExpressionTreesMod( false),
00107         m_LagrangianExpTree(NULL),
00108         m_bLagrangianExpTreeCreated( false),
00109         m_LagrangianSparseHessian( NULL),
00110         m_bLagrangianSparseHessianCreated( false),
00111         m_miNonLinearVarsReverseMap( NULL),
00112         m_bAllNonlinearVariablesIndex( false),
00113         m_bOSADFunIsCreated( false),
00114         m_bCppADTapesBuilt( false),
00115         m_bCppADMustReTape( false),
00116         m_bDuplicateExpressionTreesMap( false),
00117         m_bNonLinearStructuresInitialized( false),
00118         m_bSparseJacobianCalculated( false),
00119         m_iHighestOrderEvaluated( -1),
00120         m_mmdObjGradient( NULL),
00121         m_bProcessTimeDomain( false),
00122         m_bProcessTimeStages( false),
00123         m_bProcessTimeInterval( false),
00124         m_bFiniteTimeStages( false),
00125         m_iNumberOfTimeStages(-1),
00126         m_msTimeDomainStageNames(NULL),
00127         m_miTimeDomainStageVariableNumber(NULL),
00128         m_mmiTimeDomainStageVarList(NULL),
00129         m_miTimeDomainStageConstraintNumber(NULL),
00130         m_mmiTimeDomainStageConList(NULL),
00131         m_miTimeDomainStageObjectiveNumber(NULL),
00132         m_mmiTimeDomainStageObjList(NULL),
00133         bUseExpTreeForFunEval( false)
00134 
00135 {    
00136         #ifdef DEBUG
00137         cout << "Inside OSInstance Constructor" << endl;
00138         #endif
00139         this->instanceHeader = new InstanceHeader();
00140         this->instanceData = new InstanceData();
00141 }  
00142 
00143 OSInstance::~OSInstance(){
00144         #ifdef DEBUG
00145         cout << "OSInstance Destructor Called" << endl;
00146         #endif
00147         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
00148         // delete  the temporary arrays
00149         
00150         if(m_bProcessVariables == true){
00151                 delete[] m_msVariableNames;
00152                 m_msVariableNames = NULL;
00153                 delete[] m_mcVariableTypes;
00154                 m_mcVariableTypes = NULL;
00155                 delete[] m_mdVariableLowerBounds;
00156                 m_mdVariableLowerBounds = NULL;
00157                 delete[] m_mdVariableUpperBounds;
00158                 m_mdVariableUpperBounds = NULL;
00159         }
00160         
00161         
00162         if(m_bProcessConstraints == true){
00163                 delete[] m_msConstraintNames;
00164                 m_msConstraintNames = NULL;
00165                 delete[] m_mcConstraintTypes;
00166                 m_mcConstraintTypes = NULL;
00167                 delete[]  m_mdConstraintConstants;
00168                 m_mdConstraintConstants = NULL;
00169                 delete[] m_mdConstraintLowerBounds;
00170                 m_mdConstraintLowerBounds = NULL;
00171                 delete[] m_mdConstraintUpperBounds;
00172                 m_mdConstraintUpperBounds = NULL;
00173         }
00174         
00175         
00176         int i;
00177         //if(instanceData->objectives->numberOfObjectives > 0 && m_mObjectiveCoefficients != NULL){
00178         if(m_bProcessObjectives == true ){
00179                 for(i = 0; i < instanceData->objectives->numberOfObjectives; i++){
00180                         #ifdef DEBUG
00181                         std::cout <<  "Delete m_mObjectiveCoefficients[i]" << std::endl;
00182                         #endif
00183                         delete m_mObjectiveCoefficients[i];
00184                         m_mObjectiveCoefficients[i] = NULL;
00185                 }
00186                 #ifdef DEBUG
00187                 std::cout <<  "Delete m_msObjectiveNames" << std::endl;
00188                 std::cout <<  "Delete m_msMaxOrMins" << std::endl;
00189                 std::cout <<  "Delete m_miNumberOfObjCoef" << std::endl;
00190                 std::cout <<  "Delete m_mdObjectiveConstants" << std::endl;
00191                 std::cout <<  "Delete m_mdObjectiveWeights" << std::endl;
00192                 #endif          
00193                 delete[] m_msObjectiveNames;
00194                 m_msObjectiveNames = NULL;
00195                 delete[] m_msMaxOrMins;
00196                 m_msMaxOrMins = NULL;
00197                 delete[] m_miNumberOfObjCoef;
00198                 m_miNumberOfObjCoef = NULL;
00199                 delete[] m_mdObjectiveConstants; 
00200                 m_mdObjectiveConstants = NULL;
00201                 delete[] m_mdObjectiveWeights;
00202                 m_mdObjectiveWeights = NULL;
00203                 delete[] m_mObjectiveCoefficients;
00204                 m_mObjectiveCoefficients = NULL;
00205         }
00206         
00207         if(m_bGetDenseObjectives == true){
00208                 for(i = 0; i < instanceData->objectives->numberOfObjectives; i++){
00209                         //delete m_mmdDenseObjectiveCoefficients[i];
00210                         #ifdef DEBUG
00211                         std::cout <<  "delete m_mmdDenseObjectiveCoefficients[i]" << std::endl;
00212                         #endif
00213                     delete[] m_mmdDenseObjectiveCoefficients[i];
00214                         m_mmdDenseObjectiveCoefficients[i] = NULL;
00215                 }
00216                 delete[] m_mmdDenseObjectiveCoefficients;
00217                 m_mmdDenseObjectiveCoefficients = NULL;
00218         }
00219         
00220         
00221         
00222         //if(m_bProcessLinearConstraintCoefficients == true && m_bColumnMajor == true) delete m_linearConstraintCoefficientsInColumnMajor;
00223         //if(m_bProcessLinearConstraintCoefficients == true && m_bColumnMajor == false) delete m_linearConstraintCoefficientsInRowMajor;
00224         
00225         if( m_linearConstraintCoefficientsInColumnMajor != NULL) delete m_linearConstraintCoefficientsInColumnMajor;
00226         if (m_linearConstraintCoefficientsInRowMajor != NULL ) delete m_linearConstraintCoefficientsInRowMajor;
00227         
00228         
00229         if( (m_binitForAlgDiff == true)  ){     
00230                 delete[] m_miNonLinearVarsReverseMap;
00231                 m_miNonLinearVarsReverseMap = NULL;
00232                 if(instanceData->objectives->numberOfObjectives > 0 && m_mmdObjGradient != NULL){
00233                         
00234                         #ifdef DEBUG
00235                         std::cout <<  "The number of objectives =  " << instanceData->objectives->numberOfObjectives << std::endl;
00236                         #endif
00237                         for(i = 0; i < instanceData->objectives->numberOfObjectives; i++){
00238                                 #ifdef DEBUG
00239                                 std::cout << "deleting Objective function gradient " << i << std::endl;
00240                                 #endif
00241                                 delete[] m_mmdObjGradient[i];
00242 
00243                                 m_mmdObjGradient[i] = NULL;
00244                         }
00245                         delete[] m_mmdObjGradient;
00246                         m_mmdObjGradient = NULL;
00247                 }
00248         }
00249 
00250         //std::cout << "Do garbage collection for the nonlinear API" << std::endl;
00251         // garbage collection for the gradient
00252         if(m_bNonLinearStructuresInitialized == true ){
00253                 delete[] m_mdObjectiveFunctionValues; 
00254                 m_mdObjectiveFunctionValues = NULL;     
00255                 delete[] m_mdConstraintFunctionValues;
00256                 m_mdConstraintFunctionValues = NULL;
00257         }
00258         if(m_bSparseJacobianCalculated == true){
00259                 delete[] m_miJacStart;
00260                 m_miJacStart = NULL;
00261                 delete[] m_miJacIndex;
00262                 m_miJacIndex = NULL;
00263                 delete[] m_mdJacValue;
00264                 m_mdJacValue = NULL;
00265                 delete[] m_miJacNumConTerms;
00266                 m_miJacNumConTerms = NULL;
00267         }
00268         if( m_bLagrangianExpTreeCreated == true){
00269                 delete m_LagrangianExpTree;
00270                 m_LagrangianExpTree = NULL;
00271         }
00272         if( m_bLagrangianSparseHessianCreated == true){
00273                 delete m_LagrangianSparseHessian;
00274                 m_LagrangianSparseHessian = NULL;
00275         }
00276         if( m_bSparseJacobianCalculated == true){
00277                 delete m_sparseJacMatrix;
00278                 m_sparseJacMatrix = NULL;
00279         }
00280         if( (instanceData->quadraticCoefficients->qTerm != NULL) && (m_bProcessQuadraticTerms == true) ){
00281                 delete m_quadraticTerms;
00282                 m_quadraticTerms = NULL;
00283         }
00284         if( (instanceData->quadraticCoefficients->qTerm != NULL) && (m_bQuadraticRowIndexesProcessed == true) ){
00285                 delete[] m_miQuadRowIndexes;
00286                 m_miQuadRowIndexes = NULL;
00287         }
00288         //
00289         // delete the new expression trees that got created
00290         //if( m_bLagrangianExpTreeCreated == false  ||  m_bLagrangianExpTreeCreated == true){
00291         if( (m_bProcessExpressionTrees == true) && (m_bDuplicateExpressionTreesMap == false)  ) {
00292                 for(posMapExpTree = m_mapExpressionTrees.begin(); posMapExpTree != m_mapExpressionTrees.end(); ++posMapExpTree){
00293                         std::cout << "Deleting an expression tree from the map for row  " << posMapExpTree->first  << std::endl;
00294                         delete m_mapExpressionTrees[ posMapExpTree->first ];
00295                 }
00296         }
00297         if( m_bDuplicateExpressionTreesMap == true)   {
00298                 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){          
00299                         #ifdef DEBUG
00300                                 std::cout << "Deleting an expression tree from m_mapExpressionTreesMod" << std::endl;
00301                         #endif
00302                         delete m_mapExpressionTreesMod[ posMapExpTree->first ];
00303                 }
00304         }
00305         //}
00307         if( (m_bNonlinearExpressionTreeIndexesProcessed == true) && (m_mapExpressionTrees.size() > 0) ){
00308                 std::cout << "Deleting  m_miNonlinearExpressionTreeIndexes" << std::endl;
00309                 delete[] m_miNonlinearExpressionTreeIndexes;
00310                 std::cout << "Done Deleting  m_miNonlinearExpressionTreeIndexes" << std::endl;
00311                 m_miNonlinearExpressionTreeIndexes = NULL;
00312         }
00313         if( (m_bNonlinearExpressionTreeModIndexesProcessed == true) && (m_mapExpressionTreesMod.size() > 0) ){
00314                 std::cout << "Deleting  m_miNonlinearExpressionTreeModIndexes" << std::endl;
00315                 delete[] m_miNonlinearExpressionTreeModIndexes;
00316                 std::cout << "Done Deleting  m_miNonlinearExpressionTreeModIndexes" << std::endl;
00317                 m_miNonlinearExpressionTreeModIndexes = NULL;
00318         }
00319         if(m_bOSADFunIsCreated == true){
00320                 try{
00321 #ifdef COIN_HAS_CPPAD
00322                         delete Fad;
00323                         Fad = NULL;
00324 #else
00325                 throw ErrorClass( "Error: An Algorithmic Differentiation Package Not Available");
00326 #endif
00327                 }
00328                 catch(const ErrorClass& eclass){
00329                         throw ErrorClass( eclass.errormsg);
00330                 }
00331         }
00332 //      if( (instanceData->timeDomain->stages->stage != NULL) && (m_bProcessTimeStages == true) ){
00333 //              delete m_Stages;
00334 //              m_Stages = NULL;
00335 //      }
00336 
00337         if (m_msTimeDomainStageNames != NULL) {
00338                 delete[] m_msTimeDomainStageNames;
00339                 m_msTimeDomainStageNames = NULL;
00340         }
00341 
00342         if (m_miTimeDomainStageVariableNumber != NULL) {
00343                 delete[] m_miTimeDomainStageVariableNumber;
00344                 m_miTimeDomainStageVariableNumber = NULL;
00345         }
00346 
00347         if (m_mmiTimeDomainStageVarList != NULL) {
00348                 for (int i = 0; i < m_iNumberOfTimeStages; i ++) 
00349                         delete[] m_mmiTimeDomainStageVarList[i];
00350                 delete[] m_mmiTimeDomainStageVarList;
00351                 m_mmiTimeDomainStageVarList = NULL;
00352         }
00353 
00354         if (m_miTimeDomainStageConstraintNumber != NULL) {
00355                 delete[] m_miTimeDomainStageConstraintNumber;
00356                 m_miTimeDomainStageConstraintNumber = NULL;
00357         }
00358 
00359         if (m_mmiTimeDomainStageConList != NULL) {
00360                 for (int i = 0; i < m_iNumberOfTimeStages; i ++) 
00361                         delete[] m_mmiTimeDomainStageConList[i];
00362                 delete[] m_mmiTimeDomainStageConList;
00363                 m_mmiTimeDomainStageConList = NULL;
00364         }
00365 
00366         if (m_miTimeDomainStageObjectiveNumber != NULL) {
00367                 delete[] m_miTimeDomainStageObjectiveNumber;
00368                 m_miTimeDomainStageObjectiveNumber = NULL;
00369         }
00370 
00371         if (m_mmiTimeDomainStageObjList != NULL) {
00372                 for (int i = 0; i < m_iNumberOfTimeStages; i ++) 
00373                         delete[] m_mmiTimeDomainStageObjList[i];
00374                 delete[] m_mmiTimeDomainStageObjList;
00375                 m_mmiTimeDomainStageObjList = NULL;
00376         }
00377 
00378         // delete the two children of OSInstance
00379         //delete instanceHeader object
00380         delete instanceHeader;
00381         instanceHeader = NULL;
00382         //delete instanceData object
00383         delete instanceData;
00384         instanceData = NULL;
00385 }//OSInstance Destructor
00386 
00387 InstanceHeader::InstanceHeader():
00388         description(""),
00389         name(""),
00390         source("")
00391 
00392 { 
00393         #ifdef DEBUG
00394         cout << "Inside the InstanceHeader Constructor" << endl;
00395         #endif
00396 } 
00397 
00398 
00399 InstanceHeader::~InstanceHeader(){
00400         #ifdef DEBUG  
00401         cout << "Inside the InstanceHeader Destructor" << endl;
00402         #endif
00403 } 
00404 
00405 Variable::Variable():
00406         lb(0.0),
00407         ub(OSDBL_MAX),
00408         //init(OSNAN),  deprecated
00409         type('C'), 
00410         name("")
00411         //initString("") deprecated
00412 {  
00413         #ifdef DEBUG
00414         cout << "Inside the Variable Constructor" << endl;
00415         #endif 
00416 } 
00417 
00418 Variable::~Variable(){  
00419         #ifdef DEBUG
00420         cout << "Inside the Variable Destructor" << endl; 
00421         #endif
00422 } 
00423 
00424 Variables::Variables(){  
00425         #ifdef DEBUG
00426         cout << "Inside the Variables Constructor" << endl; 
00427         #endif 
00428         numberOfVariables = 0;
00429         var = NULL; 
00430 }
00431 
00432 Variables::~Variables(){ 
00433         #ifdef DEBUG 
00434         cout << "Inside the Variables Destructor" << endl;
00435         #endif
00436         int i;
00437         if(numberOfVariables > 0 && var != NULL){
00438                 for(i = 0; i < numberOfVariables; i++){
00439                         #ifdef DEBUG 
00440                         cout << "Deleting var[ i]" << endl;
00441                         #endif
00442                         delete var[i];
00443                         var[i] = NULL;
00444                 }
00445         }
00446         delete[] var;
00447         var = NULL; 
00448 }  
00449 
00450 ObjCoef::ObjCoef():
00451         idx(-1),
00452         value(0.0)  
00453 {  
00454         #ifdef DEBUG
00455         cout << "Inside the Coef Constructor" << endl;
00456         #endif 
00457 }
00458 
00459 ObjCoef::~ObjCoef(){ 
00460         #ifdef DEBUG
00461         cout << "Inside the ObjCoef Destructor" << endl;  
00462         #endif
00463 }
00464 
00465 Objective::Objective():
00466         name("") ,
00467         maxOrMin("min"),
00468         constant(0.0),
00469         weight(1.0),
00470         numberOfObjCoef(0),
00471         coef(NULL)
00472 { 
00473  
00474         #ifdef DEBUG
00475         cout << "Inside the Objective Constructor" << endl;
00476         #endif
00477 }
00478 
00479 Objective::~Objective(){
00480         #ifdef DEBUG  
00481         cout << "Inside the Objective Destructor" << endl;
00482         #endif
00483         int i;
00484         if(numberOfObjCoef > 0 && coef != NULL){
00485                 for(i = 0; i < numberOfObjCoef; i++){
00486                         delete coef[i];
00487                         coef[i] = NULL;
00488                 }
00489         }
00490         delete[] coef;
00491         coef = NULL;
00492 }  
00493 
00494 Objectives::Objectives()
00495 {  
00496         #ifdef DEBUG
00497         cout << "Inside the Objectives Constructor" << endl; 
00498         #endif
00499         numberOfObjectives = 0;
00500         obj = NULL;
00501 } 
00502 
00503 Objectives::~Objectives(){ 
00504         #ifdef DEBUG 
00505         cout << "Inside the Objectives Destructor" << endl;
00506         #endif
00507         int i;
00508         if(numberOfObjectives > 0 && obj != NULL){
00509                 for(i = 0; i < numberOfObjectives; i++){
00510                         delete obj[i];
00511                         obj[i] = NULL;
00512                 }
00513         }
00514         delete[] obj;
00515         obj = NULL;
00516 }
00517 
00518 Constraint::Constraint():
00519         name(""),
00520         constant(0.0),
00521         lb(-OSDBL_MAX),
00522         ub(OSDBL_MAX)
00523 
00524 {
00525         #ifdef DEBUG  
00526         cout << "Inside the Constraint Constructor" << endl;
00527         #endif
00528 } 
00529 
00530 Constraint::~Constraint(){  
00531         #ifdef DEBUG
00532         cout << "Inside the Constraint Destructor" << endl;
00533         #endif
00534 } 
00535 
00536 Constraints::Constraints():
00537         numberOfConstraints(0),
00538         con(NULL)
00539 {
00540         #ifdef DEBUG
00541         cout << "Inside the Constraints Constructor" << endl;
00542         #endif
00543 } 
00544 
00545 Constraints::~Constraints(){  
00546         #ifdef DEBUG
00547         cout << "Inside the Constraints Destructor" << endl;
00548         #endif
00549         int i;
00550         if(numberOfConstraints > 0 && con != NULL){
00551                 for( i = 0; i < numberOfConstraints; i++){
00552                         delete con[i];
00553                         con[i] = NULL;
00554                 }
00555         }
00556         delete[] con;
00557         con = NULL;
00558 } 
00559 
00560 
00561 
00562 LinearConstraintCoefficients::LinearConstraintCoefficients():
00563         numberOfValues(0) ,
00564         iNumberOfStartElements( 0)
00565 { 
00566         #ifdef DEBUG 
00567         cout << "Inside the LinearConstraintCoefficients Constructor" << endl; 
00568         #endif
00569         start = new IntVector();
00570         rowIdx = new IntVector();
00571         colIdx = new IntVector();
00572         value = new DoubleVector();
00573 
00574 } 
00575 
00576 
00577 LinearConstraintCoefficients::~LinearConstraintCoefficients(){  
00578         #ifdef DEBUG
00579         cout << "Inside the LinearConstraintCoefficients Destructor" << endl; 
00580         #endif
00581         delete start;
00582         start = NULL;
00583         delete rowIdx;
00584         rowIdx = NULL;
00585         delete colIdx;
00586         colIdx = NULL;
00587         delete value;
00588         value = NULL;
00589 }
00590 
00591 QuadraticTerm::QuadraticTerm():
00592 
00593         idx(0),   
00594         idxOne(-1),
00595         idxTwo(-1),
00596         coef(0.0)
00597 
00598 {
00599         #ifdef DEBUG  
00600         cout << "Inside the QuadraticTerm Constructor" << endl;
00601         #endif
00602 } 
00603 
00604 
00605 QuadraticTerm::~QuadraticTerm(){  
00606         #ifdef DEBUG
00607         cout << "Inside the QuadraticTerm Destructor" << endl;
00608         #endif
00609 }
00610 
00611 
00612 
00613 QuadraticCoefficients::QuadraticCoefficients():
00614         numberOfQuadraticTerms(0),
00615         qTerm(NULL)
00616 { 
00617         #ifdef DEBUG 
00618         cout << "Inside the QuadraticCoefficients Constructor" << endl;
00619         #endif
00620 }//end QuadraticCoefficients() 
00621 
00622 
00623 QuadraticCoefficients::~QuadraticCoefficients(){
00624         #ifdef DEBUG  
00625         cout << "Inside the QuadraticCoefficients Destructor" << endl;
00626         #endif
00627         int i;
00628         if(numberOfQuadraticTerms > 0 && qTerm != NULL){
00629                 for( i = 0; i < numberOfQuadraticTerms; i++){
00630                         delete qTerm[i];
00631                         qTerm[i] = NULL;
00632                 }
00633         }
00634         delete[] qTerm;
00635         qTerm = NULL;  
00636 }//end ~QuadraticCoefficients()  
00637 
00638 
00639 Nl::Nl(){
00640         idx = 0;
00641         osExpressionTree = NULL;
00642         m_bDeleteExpressionTree = true;
00643 }//end Nl
00644  
00645  
00646 Nl::~Nl(){
00647         #ifdef DEBUG  
00648         cout << "Inside the Nl Destructor" << endl;
00649         #endif
00650         // don't delete the expression tree if we created a map of the expression
00651         // trees, otherwise we would destroy twice
00652         if( m_bDeleteExpressionTree == true){
00653                 delete osExpressionTree;
00654                 osExpressionTree = NULL;
00655         }
00656 }//end ~Nl
00657 
00658 
00659 
00660 NonlinearExpressions::NonlinearExpressions():
00661         numberOfNonlinearExpressions(0) ,
00662         nl(NULL)
00663 { 
00664         #ifdef DEBUG 
00665         cout << "Inside the NonlinearExpressions Constructor" << endl;
00666         #endif
00667 }//end NonlinearExpressions() 
00668 
00669 NonlinearExpressions::~NonlinearExpressions(){
00670         #ifdef DEBUG  
00671         cout << "Inside the NonlinearExpressions Destructor" << endl;
00672         cout << "NUMBER OF NONLINEAR EXPRESSIONS = " << numberOfNonlinearExpressions << endl;
00673         #endif
00674         int i;
00675         if(numberOfNonlinearExpressions > 0 && nl != NULL){
00676                 for( i = 0; i < numberOfNonlinearExpressions; i++){
00677                         #ifdef DEBUG  
00678                                 cout << "DESTROYING EXPRESSION " << nl[ i]->idx << endl;
00679                         #endif
00680                         delete nl[i];
00681                         nl[i] = NULL;
00682                 }
00683         }
00684         delete[] nl;
00685         nl = NULL;  
00686 }//end ~NonlinearExpressions()  
00687 
00688 
00689 TimeDomainStageVar::TimeDomainStageVar():
00690         idx(0)
00691 {
00692         #ifdef DEBUG 
00693         cout << "Inside the Stage Objectives Var Constructor" << endl;
00694         #endif
00695 } // end TimeDomainStageVar
00696 
00697 
00698 TimeDomainStageVar::~TimeDomainStageVar()
00699 {
00700         #ifdef DEBUG 
00701         cout << "Inside the Stage Objectives Var Destructor" << endl;
00702         #endif
00703 } // end ~TimeDomainStageVar
00704 
00705 
00706 TimeDomainStageVariables::TimeDomainStageVariables():
00707         numberOfVariables(0),
00708         startIdx(-1)
00709 {
00710         #ifdef DEBUG 
00711         cout << "Inside the Stage Variables Constructor" << endl;
00712         #endif
00713         var = NULL;
00714 } // end TimeDomainStageVariables
00715 
00716 TimeDomainStageVariables::~TimeDomainStageVariables()
00717 {
00718         #ifdef DEBUG 
00719         cout << "Inside the Stage Variables Destructor" << endl;
00720         #endif
00721         if (numberOfVariables > 0 && var != NULL){
00722                 for (int i = 0; i < numberOfVariables; i++) {
00723                         delete var[i];
00724                         var[i] = NULL;
00725                 }
00726         }
00727         delete [] var;
00728         var = NULL;
00729 } // end ~TimeDomainStageVariables
00730 
00731 TimeDomainStageCon::TimeDomainStageCon():
00732         idx(0)
00733 {
00734         #ifdef DEBUG 
00735         cout << "Inside the Stage Objectives Con Constructor" << endl;
00736         #endif
00737 } // end TimeDomainStageCon
00738 
00739 
00740 TimeDomainStageCon::~TimeDomainStageCon()
00741 {
00742         #ifdef DEBUG 
00743         cout << "Inside the Stage Objectives Con Destructor" << endl;
00744         #endif
00745 } // end ~TimeDomainStageCon
00746 
00747 
00748 TimeDomainStageConstraints::TimeDomainStageConstraints():
00749         numberOfConstraints(0),
00750         startIdx(-1)
00751 {
00752         #ifdef DEBUG 
00753         cout << "Inside the Stage Constraints Constructor" << endl;
00754         #endif
00755         con = NULL;
00756 } // end TimeDomainStageConstraints
00757 
00758 TimeDomainStageConstraints::~TimeDomainStageConstraints()
00759 {
00760         #ifdef DEBUG 
00761         cout << "Inside the Stage Constraints Destructor" << endl;
00762         #endif
00763         if (numberOfConstraints > 0 && con != NULL){
00764                 for (int i = 0; i < numberOfConstraints; i++) {
00765                         delete con[i];
00766                         con[i] = NULL;
00767                 }
00768         }
00769         delete [] con;
00770         con = NULL;
00771 } // end ~TimeDomainStageConstraints
00772 
00773 TimeDomainStageObj::TimeDomainStageObj():
00774         idx(0)
00775 {
00776         #ifdef DEBUG 
00777         cout << "Inside the Stage Objectives Obj Constructor" << endl;
00778         #endif
00779 } // end TimeDomainStageObj
00780 
00781 
00782 TimeDomainStageObj::~TimeDomainStageObj()
00783 {
00784         #ifdef DEBUG 
00785         cout << "Inside the Stage Objectives Obj Destructor" << endl;
00786         #endif
00787 } // end ~TimeDomainStageObj
00788 
00789 
00790 TimeDomainStageObjectives::TimeDomainStageObjectives():
00791         numberOfObjectives(0),
00792         startIdx(-1)
00793 {
00794         #ifdef DEBUG 
00795         cout << "Inside the Stage Objectives Constructor" << endl;
00796         #endif
00797         obj = NULL;
00798 } // end TimeDomainStageObjectives
00799 
00800 TimeDomainStageObjectives::~TimeDomainStageObjectives()
00801 {
00802         #ifdef DEBUG 
00803         cout << "Inside the Stage Objectives Destructor" << endl;
00804         #endif
00805         if (numberOfObjectives > 0 && obj != NULL){
00806                 for (int i = 0; i < numberOfObjectives; i++) {
00807                         delete obj[i];
00808                         obj[i] = NULL;
00809                 }
00810         }
00811         delete [] obj;
00812         obj = NULL;
00813 } // end ~TimeDomainStageObjectives
00814 
00815 TimeDomainStage::TimeDomainStage():
00816         name("")
00817 { 
00818         #ifdef DEBUG 
00819         cout << "Inside the Stage Constructor" << endl;
00820         #endif
00821         variables   = NULL;
00822         constraints = NULL;
00823         objectives  = NULL;
00824 }//end TimeDomainStage() 
00825 
00826 
00827 TimeDomainStage::~TimeDomainStage(){
00828         #ifdef DEBUG  
00829         cout << "Inside the Stage Destructor" << endl;
00830         #endif
00831         if (variables != NULL)
00832         {       delete variables;
00833                 variables = NULL;
00834         }
00835         if (constraints != NULL)
00836         {       delete constraints;
00837                 constraints = NULL;
00838         }
00839         if (objectives != NULL)
00840         {       delete objectives;
00841                 objectives = NULL;
00842         }
00843 }//end ~TimeDomainStage()  
00844 
00845 
00846 TimeDomainStages::TimeDomainStages():
00847         numberOfStages(0),
00848         stage(NULL)
00849 {
00850         #ifdef DEBUG  
00851         cout << "Inside the Stages Constructor" << endl;
00852         #endif
00853 } 
00854 
00855 
00856 TimeDomainStages::~TimeDomainStages(){  
00857         #ifdef DEBUG
00858         cout << "Inside the Stages Destructor" << endl;
00859         #endif
00860         int i;
00861         if(numberOfStages > 0 && stage != NULL){
00862                 for( i = 0; i < numberOfStages; i++){
00863                         delete stage[i];
00864                         stage[i] = NULL;
00865                 }
00866         }
00867         delete[] stage;
00868         stage = NULL;  
00869 }
00870 
00871 TimeDomainInterval::TimeDomainInterval():
00872         start(0.0),
00873         horizon(0.0)
00874 {
00875         #ifdef DEBUG  
00876         cout << "Inside the Interval Constructor" << endl;
00877         #endif
00878 } 
00879 
00880 
00881 TimeDomainInterval::~TimeDomainInterval(){  
00882         #ifdef DEBUG
00883         cout << "Inside the Interval Destructor" << endl;
00884         #endif
00885 }
00886 
00887 TimeDomain::TimeDomain()
00888 {
00889         #ifdef DEBUG
00890         cout << "Inside the TimeDomain Constructor" << endl;
00891         #endif
00892         stages = NULL;
00893         interval = NULL;
00894 }
00895 
00896 TimeDomain::~TimeDomain()
00897 {  
00898         #ifdef DEBUG
00899         cout << "Inside the TimeDomain Destructor" << endl;
00900         #endif
00901         if (stages != NULL)
00902         {       delete stages;
00903                 stages = NULL;
00904         }
00905         if (interval != NULL)
00906         {       delete interval;
00907                 interval = NULL;
00908         }
00909 } 
00910 
00911 
00912 InstanceData::InstanceData(){ 
00913         #ifdef DEBUG 
00914         cout << "Inside the InstanceData Constructor" << endl;
00915         #endif 
00916         variables = new Variables();
00917         objectives = new Objectives();
00918         constraints = new Constraints();
00919         linearConstraintCoefficients = new LinearConstraintCoefficients();
00920         quadraticCoefficients = new QuadraticCoefficients();
00921         nonlinearExpressions = new NonlinearExpressions();
00922         timeDomain = NULL;
00923 } 
00924 
00925 InstanceData::~InstanceData(){  
00926         #ifdef DEBUG
00927         cout << "Inside the InstanceData Destructor" << endl; 
00928         #endif
00929         delete variables;
00930         variables = NULL;
00931         delete objectives;
00932         objectives = NULL;
00933         delete constraints;
00934         constraints = NULL;
00935         delete linearConstraintCoefficients;
00936         linearConstraintCoefficients = NULL;
00937         delete quadraticCoefficients;
00938         quadraticCoefficients = NULL;
00939         delete nonlinearExpressions;
00940         nonlinearExpressions = NULL;
00941         if (timeDomain != NULL)
00942         {   delete timeDomain;
00943                 timeDomain = NULL;
00944         }
00945 } 
00946 
00947 string OSInstance::getInstanceName(){
00948         if(  m_sInstanceName.length() <= 0){
00949                 m_sInstanceName = instanceHeader->name;
00950         }
00951         return m_sInstanceName;
00952 }//getInstanceName
00953 
00954 
00955 string OSInstance::getInstanceSource(){
00956         if( m_sInstanceSource.length() <= 0){
00957                 m_sInstanceSource = instanceHeader->source;
00958         }
00959         return m_sInstanceSource;
00960 }//getInstanceSource
00961 
00962 string OSInstance::getInstanceDescription(){
00963         if(m_sInstanceDescription.length() <= 0){
00964                 m_sInstanceDescription = instanceHeader->description;
00965         }
00966         return m_sInstanceDescription;
00967 }//getInstanceDescription
00968 
00969 int OSInstance::getVariableNumber(){
00970         if(m_iVariableNumber == -1){
00971                 m_iVariableNumber = instanceData->variables->numberOfVariables;
00972         }
00973         return m_iVariableNumber;
00974 }//getVariableNumber
00975 
00976 int OSInstance::getNumberOfNonlinearExpressions(){
00977         if(m_iNonlinearExpressionNumber == -1){
00978                 m_iNonlinearExpressionNumber = instanceData->nonlinearExpressions->numberOfNonlinearExpressions;
00979         }
00980         return m_iNonlinearExpressionNumber;
00981 }//getNumberOfNonlinearExpressions
00982 
00983 
00984 
00985 bool OSInstance::processVariables() {
00986         if(m_bProcessVariables == true && bVariablesModified == false) return true;
00987         //m_bProcessVariables = true;
00988         string vartype ="CBIS";
00989         int i = 0;
00990         int n = instanceData->variables->numberOfVariables;
00991         try{
00992                 if(instanceData->variables->var == NULL) throw ErrorClass("no variables defined");
00993                 /*
00994                         m_mdVariableInitialValues = new double[n];
00995                         for(i = 0; i < n; i++){
00996                                 if(CommonUtil::ISOSNAN(instanceData->variables->var[ 0]->init) == true ){                               
00997                                         m_mdVariableInitialValues[i] =  1.7171;
00998                                 }
00999                                 else{
01000                                         m_mdVariableInitialValues[i] = instanceData->variables->var[i]->init;
01001                                 }
01002                         }
01003                 //}
01004                 if((instanceData->variables->var[0]->initString.length() > 0)){
01005                         m_msVariableInitialStringValues = new string[n];
01006                         for(i = 0; i < n; i++) m_msVariableInitialStringValues[i] = instanceData->variables->var[i]->initString;
01007                 }
01008                 */
01009                 
01010                 if(m_bProcessVariables != true ){
01011                         m_mcVariableTypes = new char[n];
01012                         m_mdVariableLowerBounds = new double[n];
01013                         m_mdVariableUpperBounds = new double[n];
01014                         m_msVariableNames = new string[n];
01015                         m_bProcessVariables = true;
01016                 }
01017                 m_iNumberOfBinaryVariables = 0;
01018                 m_iNumberOfIntegerVariables = 0;
01019                 m_iNumberOfStringVariables = 0;
01020                 for(i = 0; i < n; i++){
01021                         if(vartype.find(instanceData->variables->var[i]->type) == string::npos) throw ErrorClass("wrong variable type");
01022                         m_mcVariableTypes[i] = instanceData->variables->var[i]->type;
01023                         if(m_mcVariableTypes[i] == 'B') m_iNumberOfBinaryVariables++;
01024                         if(m_mcVariableTypes[i] == 'I') m_iNumberOfIntegerVariables++;
01025                         if(m_mcVariableTypes[i] == 'S') m_iNumberOfStringVariables++;
01026                         m_mdVariableLowerBounds[i] = instanceData->variables->var[i]->lb;
01027                         m_mdVariableUpperBounds[i] = instanceData->variables->var[i]->ub;
01028                 }
01029                 if(instanceData->variables->var[0]->name.length() > 0 || instanceData->variables->var[n-1]->name.length() > 0){
01030                         for(i = 0; i < n; i++) m_msVariableNames[i] = instanceData->variables->var[i]->name;
01031                 } 
01032                 return true;
01033         } //end try
01034         catch(const ErrorClass& eclass){
01035                 throw ErrorClass( eclass.errormsg);
01036         } 
01037 }//processVariables
01038         
01039 string* OSInstance::getVariableNames() {
01040         processVariables();
01041         return m_msVariableNames;
01042 }//getVariableNames     
01043 
01044 /*
01045 double* OSInstance::getVariableInitialValues() {
01046         processVariables();
01047         return m_mdVariableInitialValues;
01048 }//getVariableInitialValues
01049 
01050 string* OSInstance::getVariableInitialStringValues() {
01051         processVariables();
01052         return m_msVariableInitialStringValues;
01053 }//getVariableInitialStringValues
01054 */
01055 
01056 char* OSInstance::getVariableTypes() {
01057         processVariables();
01058         return m_mcVariableTypes;
01059 }//getVariableTypes
01060 
01061 int OSInstance::getNumberOfIntegerVariables() {
01062         processVariables();
01063         return m_iNumberOfIntegerVariables;
01064 }//getNumberOfIntegerVariables
01065 
01066 int OSInstance::getNumberOfBinaryVariables() {
01067         processVariables();
01068         return m_iNumberOfBinaryVariables;
01069 }//getNumberOfBinaryVariables
01070 
01071 int OSInstance::getNumberOfStringVariables() {
01072         processVariables();
01073         return m_iNumberOfStringVariables;
01074 }//getNumberOfStringVariables
01075 
01076 double* OSInstance::getVariableLowerBounds() {
01077         processVariables();
01078         return m_mdVariableLowerBounds;
01079 }//getVariableLowerBounds
01080 
01081 double* OSInstance::getVariableUpperBounds() {
01082         processVariables();
01083         return m_mdVariableUpperBounds;
01084 }//getVariableUpperBounds
01085 
01086 int OSInstance::getObjectiveNumber(){
01087         if(m_iObjectiveNumber == -1){
01088                 m_iObjectiveNumber = instanceData->objectives->numberOfObjectives;
01089         }
01090         return m_iObjectiveNumber;
01091 }//getObjectiveNumber
01092 
01093 
01094 bool OSInstance::processObjectives() {
01095         if(m_bProcessObjectives == true && bObjectivesModified == false) return true;
01096         //m_bProcessObjectives = true;
01097         int i = 0;
01098         int j = 0;
01099         if(instanceData == NULL || instanceData->objectives == NULL || instanceData->objectives->obj == NULL || instanceData->objectives->numberOfObjectives == 0) return true;
01100         int n = instanceData->objectives->numberOfObjectives;
01101         try{
01102                 if(m_bProcessObjectives != true){
01103                         m_msMaxOrMins = new string[n];
01104                         m_miNumberOfObjCoef = new int[n];
01105                         m_mdObjectiveConstants = new double[n];
01106                         m_mdObjectiveWeights = new double[n];
01107                         m_mObjectiveCoefficients = new SparseVector*[n];
01108                         m_msObjectiveNames = new string[n];
01109                         for(i = 0; i < n; i++){
01110                                 m_mObjectiveCoefficients[i] = new SparseVector(instanceData->objectives->obj[ j]->numberOfObjCoef);
01111                                 //m_mObjectiveCoefficients[i]->bDeleteArrays=false;
01112                         }
01113                         m_bProcessObjectives = true;
01114                 }
01115 
01116                 for(i = 0; i < n; i++){
01117                         if((instanceData->objectives->obj[i]->maxOrMin.compare("max") != 0) && (instanceData->objectives->obj[i]->maxOrMin.compare("min") != 0 )) throw ErrorClass("wrong objective maxOrMin");
01118                         m_msMaxOrMins[i] = instanceData->objectives->obj[i]->maxOrMin;
01119                         m_miNumberOfObjCoef[i] = instanceData->objectives->obj[i]->numberOfObjCoef;
01120                         m_mdObjectiveConstants[i] = instanceData->objectives->obj[i]->constant;
01121                         m_mdObjectiveWeights[i] = instanceData->objectives->obj[i]->weight;
01122                         if(instanceData->objectives->obj[i]->coef == NULL && m_miNumberOfObjCoef[i] != 0){
01123                                 throw ErrorClass("objective coefficient number inconsistent with objective coefficient array");
01124                         }
01125                         if(instanceData->objectives->obj[i]->coef != NULL){
01126                                 for(j = 0; j < m_mObjectiveCoefficients[i]->number; j++){
01127                                         m_mObjectiveCoefficients[i]->indexes[j] = instanceData->objectives->obj[i]->coef[j]->idx;
01128                                         m_mObjectiveCoefficients[i]->values[j] = instanceData->objectives->obj[i]->coef[j]->value;                      
01129                                 }
01130                         }
01131                 }       
01132                 if(instanceData->objectives->obj[0]->name.length() > 0 || instanceData->objectives->obj[n-1]->name.length() > 0){                       
01133                         for(i = 0; i < n; i++) m_msObjectiveNames[i] = instanceData->objectives->obj[i]->name;
01134                 }       
01135                 return true;
01136         }
01137         catch(const ErrorClass& eclass){
01138                 throw ErrorClass( eclass.errormsg);
01139         }
01140 }//processObjectives
01141 
01142 string* OSInstance::getObjectiveNames() {
01143         processObjectives();
01144         return m_msObjectiveNames;
01145 }//getObjectiveNames
01146 
01147 string* OSInstance::getObjectiveMaxOrMins() {
01148         processObjectives();
01149         return m_msMaxOrMins;
01150 }//getObjectiveMaxOrMins
01151 
01152 int* OSInstance::getObjectiveCoefficientNumbers(){
01153                 processObjectives();
01154                 return m_miNumberOfObjCoef;
01155         }//getObjectiveCoefficientNumbers
01156 
01157 double* OSInstance::getObjectiveConstants() {
01158         processObjectives();
01159         return m_mdObjectiveConstants;
01160 }//getObjectiveConstants
01161 
01162 double* OSInstance::getObjectiveWeights() {
01163         processObjectives();
01164         return m_mdObjectiveWeights;
01165 }//getObjectiveWeights
01166 
01167 SparseVector** OSInstance::getObjectiveCoefficients() {
01168         processObjectives();
01169         return m_mObjectiveCoefficients;
01170 }//getObjectiveCoefficients
01171 
01172 
01173 double** OSInstance::getDenseObjectiveCoefficients() {
01174         if(m_bGetDenseObjectives == true && bObjectivesModified == false) return m_mmdDenseObjectiveCoefficients;
01175         //m_bGetDenseObjectives = true;
01176         int i, j, numobjcoef;
01177         SparseVector *sparsevec;
01178         if(instanceData->objectives->obj == NULL || instanceData->objectives->numberOfObjectives == 0) return m_mmdDenseObjectiveCoefficients;
01179         int m = instanceData->objectives->numberOfObjectives;
01180         int n = instanceData->variables->numberOfVariables;
01181         if(m_bGetDenseObjectives != true){
01182                 m_mmdDenseObjectiveCoefficients = new double*[m];
01183                 for(i = 0; i < m; i++){
01184                         m_mmdDenseObjectiveCoefficients[ i] = new double[n];
01185                 }
01186                 m_bGetDenseObjectives = true;
01187         }
01188 
01189         for(i = 0; i < m; i++){
01190                 sparsevec = this->getObjectiveCoefficients()[i];
01191                 for(j = 0; j < n; j++){
01192                         m_mmdDenseObjectiveCoefficients[ i][j] = 0.0;
01193                 }
01194                 sparsevec =  this->getObjectiveCoefficients()[i];
01195                 numobjcoef = sparsevec->number;
01196                 for(j = 0; j < numobjcoef; j++){
01197                         m_mmdDenseObjectiveCoefficients[i][ sparsevec->indexes[ j]]
01198                         = sparsevec->values[ j];
01199                 }
01200         }
01201         return m_mmdDenseObjectiveCoefficients;
01202 }//getDenseObjectiveCoefficients
01203 
01204 int OSInstance::getConstraintNumber(){
01205         if(m_iConstraintNumber == -1){
01206                 m_iConstraintNumber = instanceData->constraints->numberOfConstraints;
01207         }
01208         return m_iConstraintNumber;
01209 }//getConstraintNumber
01210 
01211 bool OSInstance::processConstraints() {
01212         if(m_bProcessConstraints == true && bConstraintsModified == false) return true;
01213         //m_bProcessConstraints = true;
01214         int i = 0;
01215         ostringstream outStr;
01216         if(instanceData == NULL || instanceData->constraints == NULL || instanceData->constraints->con == NULL || instanceData->constraints->numberOfConstraints == 0) return true;
01217         int n = instanceData->constraints->numberOfConstraints;
01218         try{
01219                 if(m_bProcessConstraints != true){
01220                         m_mdConstraintLowerBounds = new double[n];
01221                         m_mdConstraintUpperBounds = new double[n];
01222                         m_mdConstraintConstants = new double[n];
01223                         m_mcConstraintTypes = new char[n];
01224                         m_msConstraintNames = new string[n];
01225                         m_bProcessConstraints = true;
01226                 }
01227                 for(i = 0; i < n; i++){
01228                         m_mdConstraintLowerBounds[i] = instanceData->constraints->con[i]->lb;
01229                         m_mdConstraintUpperBounds[i] = instanceData->constraints->con[i]->ub;
01230                         m_mdConstraintConstants[i] = instanceData->constraints->con[i]->constant;
01231                         if(m_mdConstraintLowerBounds[i] == OSDBL_MAX || m_mdConstraintUpperBounds[i] == -OSDBL_MAX) {
01232                                 throw ErrorClass( outStr.str() );
01233                         }
01234                         else if(m_mdConstraintLowerBounds[i] > m_mdConstraintUpperBounds[i]) {
01235                                 outStr << "Constraint  " ;
01236                                 outStr << i;
01237                                 outStr << " is infeasible";
01238                                 throw ErrorClass( outStr.str());
01239                         }
01240                         else if(m_mdConstraintLowerBounds[i] == -OSDBL_MAX && m_mdConstraintUpperBounds[i] == OSDBL_MAX)
01241                                 m_mcConstraintTypes[i] = 'U';
01242                         else if(m_mdConstraintLowerBounds[i] == m_mdConstraintUpperBounds[i]) 
01243                                 m_mcConstraintTypes[i] = 'E';
01244                         else if(m_mdConstraintLowerBounds[i] == -OSDBL_MAX)
01245                                 m_mcConstraintTypes[i] = 'L';
01246                         else if(m_mdConstraintUpperBounds[i] == OSDBL_MAX)
01247                                 m_mcConstraintTypes[i] = 'G';
01248                         else m_mcConstraintTypes[i] = 'R';
01249                 }
01250                 if(instanceData->constraints->con[0]->name.length() > 0 || instanceData->constraints->con[n-1]->name.length() > 0){
01251                         for(i = 0; i < n; i++) m_msConstraintNames[i] = instanceData->constraints->con[i]->name;
01252                 }
01253                 return true;
01254         }
01255         catch(const ErrorClass& eclass){
01256                 throw ErrorClass( eclass.errormsg);
01257         }
01258 }//processConstraints
01259 
01260 
01261 string* OSInstance::getConstraintNames() {
01262         processConstraints();
01263         return m_msConstraintNames;
01264 }//getConstraintNames
01265 
01266 
01267 double* OSInstance::getConstraintLowerBounds() {
01268         processConstraints();
01269         return m_mdConstraintLowerBounds;
01270 }//getConstraintLowerBounds
01271 
01272 char* OSInstance::getConstraintTypes() {
01273         processConstraints();
01274         return m_mcConstraintTypes;
01275 }//getConstraintTypes
01276 
01277 double* OSInstance::getConstraintUpperBounds() {
01278         processConstraints();
01279         return m_mdConstraintUpperBounds;
01280 }//getConstraintUpperBounds
01281 
01282 int OSInstance::getLinearConstraintCoefficientNumber(){
01283         if(m_iLinearConstraintCoefficientNumber == -1){
01284                 m_iLinearConstraintCoefficientNumber = instanceData->linearConstraintCoefficients->numberOfValues;
01285         }
01286         return m_iLinearConstraintCoefficientNumber; 
01287 }//getLinearConstraintCoefficientNumber
01288 
01289 bool OSInstance::processLinearConstraintCoefficients() {
01290         if(m_bProcessLinearConstraintCoefficients == true && bAMatrixModified == false) return true;
01291         //m_bProcessLinearConstraintCoefficients = true;
01292         try{
01293                 int n = instanceData->linearConstraintCoefficients->numberOfValues;
01294                 //value array
01295                 if((instanceData->linearConstraintCoefficients->value == NULL ) || (n == 0) ) return true;
01296                 //index array
01297                 if((instanceData->linearConstraintCoefficients->colIdx != NULL && instanceData->linearConstraintCoefficients->colIdx->el != NULL) 
01298                 && (instanceData->linearConstraintCoefficients->rowIdx != NULL && instanceData->linearConstraintCoefficients->rowIdx->el != NULL))
01299                         throw ErrorClass("ambiguous linear constraint coefficient major");
01300                 else if(instanceData->linearConstraintCoefficients->value->el == NULL) return true;
01301                 else{
01302                         if(instanceData->linearConstraintCoefficients->rowIdx->el != NULL){
01303                                 m_bColumnMajor = true;
01304                                 if(m_bProcessLinearConstraintCoefficients != true){
01305                                         m_linearConstraintCoefficientsInColumnMajor = new SparseMatrix();
01306                                         m_linearConstraintCoefficientsInColumnMajor->bDeleteArrays = false;
01307                                         m_bProcessLinearConstraintCoefficients = true;
01308                                 }
01309                                 m_linearConstraintCoefficientsInColumnMajor->isColumnMajor = true;
01310                                 m_linearConstraintCoefficientsInColumnMajor->valueSize = n;
01311                                 m_linearConstraintCoefficientsInColumnMajor->startSize = instanceData->variables->numberOfVariables + 1;
01312                         }
01313                         else{ 
01314                                 m_bColumnMajor = false; 
01315                                 if(m_bProcessLinearConstraintCoefficients != true){
01316                                         m_linearConstraintCoefficientsInRowMajor = new SparseMatrix();
01317                                         m_linearConstraintCoefficientsInRowMajor->bDeleteArrays = false;
01318                                         m_bProcessLinearConstraintCoefficients = true;
01319                                 }
01320                                 m_linearConstraintCoefficientsInRowMajor->isColumnMajor = false;
01321                                 m_linearConstraintCoefficientsInRowMajor->valueSize = n;
01322                                 m_linearConstraintCoefficientsInRowMajor->startSize = instanceData->constraints->numberOfConstraints + 1;
01323                         }
01324                 }                       
01325                 if(m_bColumnMajor == true){
01326                         m_linearConstraintCoefficientsInColumnMajor->values = instanceData->linearConstraintCoefficients->value->el;
01327                         m_linearConstraintCoefficientsInColumnMajor->indexes = instanceData->linearConstraintCoefficients->rowIdx->el;
01328                         m_linearConstraintCoefficientsInColumnMajor->starts = instanceData->linearConstraintCoefficients->start->el;                    
01329                 }
01330                 else{
01331                         m_linearConstraintCoefficientsInRowMajor->values = instanceData->linearConstraintCoefficients->value->el;
01332                         m_linearConstraintCoefficientsInRowMajor->indexes = instanceData->linearConstraintCoefficients->colIdx->el;
01333                         m_linearConstraintCoefficientsInRowMajor->starts = instanceData->linearConstraintCoefficients->start->el;                                               
01334                 }
01335                 return true;
01336         }
01337         catch(const ErrorClass& eclass){
01338                 throw ErrorClass( eclass.errormsg);
01339         }
01340 }//processLinearConstraintCoefficients
01341 
01342 bool OSInstance::getLinearConstraintCoefficientMajor() {
01343         processLinearConstraintCoefficients();  
01344         return m_bColumnMajor;          
01345 }//getLinearConstraintCoefficientMajor
01346 
01347 SparseMatrix* OSInstance::getLinearConstraintCoefficientsInColumnMajor() {
01348         processLinearConstraintCoefficients();
01349         if(m_linearConstraintCoefficientsInColumnMajor != NULL) return m_linearConstraintCoefficientsInColumnMajor;
01350         if(!m_bColumnMajor){
01351                 if(m_linearConstraintCoefficientsInRowMajor == NULL) return NULL;
01352                 m_linearConstraintCoefficientsInColumnMajor = 
01353                         MathUtil::convertLinearConstraintCoefficientMatrixToTheOtherMajor(false,
01354                                           m_linearConstraintCoefficientsInRowMajor->startSize,
01355                                           m_linearConstraintCoefficientsInRowMajor->valueSize,
01356                                           m_linearConstraintCoefficientsInRowMajor->starts,
01357                                           m_linearConstraintCoefficientsInRowMajor->indexes,
01358                                           m_linearConstraintCoefficientsInRowMajor->values,
01359                                           getVariableNumber());
01360         }
01361         return m_linearConstraintCoefficientsInColumnMajor;             
01362 }//getLinearConstraintCoefficientsInColumnMajor
01363 
01364 SparseMatrix* OSInstance::getLinearConstraintCoefficientsInRowMajor() {
01365         processLinearConstraintCoefficients();
01366         if(m_linearConstraintCoefficientsInRowMajor != NULL) return m_linearConstraintCoefficientsInRowMajor;
01367         if(m_bColumnMajor){
01368                 if(m_linearConstraintCoefficientsInColumnMajor == NULL) return NULL;
01369                 m_linearConstraintCoefficientsInRowMajor = 
01370                 MathUtil::convertLinearConstraintCoefficientMatrixToTheOtherMajor(true,
01371                                   m_linearConstraintCoefficientsInColumnMajor->startSize,
01372                                   m_linearConstraintCoefficientsInColumnMajor->valueSize,
01373                                   m_linearConstraintCoefficientsInColumnMajor->starts,
01374                                   m_linearConstraintCoefficientsInColumnMajor->indexes,
01375                                   m_linearConstraintCoefficientsInColumnMajor->values,
01376                                   getConstraintNumber());
01377         }
01378         return m_linearConstraintCoefficientsInRowMajor; 
01379 }//getLinearConstraintCoefficientsInRowMajor
01380 
01381 
01382 int OSInstance::getNumberOfQuadraticTerms(){
01383         if(m_iQuadraticTermNumber == -1){
01384         // if m_iQuadraticTermNumber == -1 then the parser did not find any q terms so 
01385         // must new the object
01386                 if(instanceData->quadraticCoefficients == NULL)instanceData->quadraticCoefficients = new QuadraticCoefficients();
01387                 m_iQuadraticTermNumber = instanceData->quadraticCoefficients->numberOfQuadraticTerms;
01388         }
01389         return m_iQuadraticTermNumber;
01390 }//getNumberOfQuadraticTerms
01391 
01392 QuadraticTerms* OSInstance::getQuadraticTerms() {
01393         if(m_bProcessQuadraticTerms) return m_quadraticTerms;
01394         m_bProcessQuadraticTerms = true;
01395         if(instanceData->quadraticCoefficients->qTerm == 0) return 0;
01396         try{
01397                 int i = 0;
01398                 QuadraticCoefficients* quadraticCoefs = instanceData->quadraticCoefficients;
01399                 int n = quadraticCoefs->numberOfQuadraticTerms;
01400                 if(!quadraticCoefs->qTerm  && n != 0) 
01401                         throw ErrorClass("quadratic term number inconsistent with quadratic term array");               
01402                 m_quadraticTerms = new QuadraticTerms();
01403                 m_quadraticTerms->rowIndexes = new int[n];
01404                 m_quadraticTerms->varOneIndexes = new int[n];
01405                 m_quadraticTerms->varTwoIndexes = new int[n];
01406                 m_quadraticTerms->coefficients = new double[n];
01407                 for(i = 0; i < n; i++){
01408                         m_quadraticTerms->rowIndexes[i] = quadraticCoefs->qTerm[i]->idx;
01409                         m_quadraticTerms->varOneIndexes[i] = quadraticCoefs->qTerm[i]->idxOne;
01410                         m_quadraticTerms->varTwoIndexes[i] = quadraticCoefs->qTerm[i]->idxTwo;
01411                         m_quadraticTerms->coefficients[i] = quadraticCoefs->qTerm[i]->coef;
01412                 } 
01413                 return m_quadraticTerms;
01414         }
01415         catch(const ErrorClass& eclass){
01416                 throw ErrorClass( eclass.errormsg);
01417         } 
01418 }//getQuadraticTerms
01419 
01420 
01421 int* OSInstance::getQuadraticRowIndexes() {
01422         if(m_bQuadraticRowIndexesProcessed == true) return m_miQuadRowIndexes;
01423         m_bQuadraticRowIndexesProcessed = true;
01424         int n = getNumberOfQuadraticTerms();    
01425         if(n <= 0) return NULL;
01426         QuadraticTerms *qTerms = NULL;
01427         qTerms = getQuadraticTerms();
01428         std::map<int, int> foundIdx;
01429         std::map<int, int>::iterator pos;
01430         int i;
01431         try{
01432                 for(i = 0; i < n; i++){
01433                         // add the terms
01434                         foundIdx[ qTerms->rowIndexes[ i] ];      
01435                 }
01436                 // now put the term into an array
01437                 m_iNumberOfQuadraticRowIndexes = foundIdx.size();
01438                 m_miQuadRowIndexes = new int[ m_iNumberOfQuadraticRowIndexes ]  ;
01439                 i = 0;
01440                 for(pos = foundIdx.begin(); pos != foundIdx.end(); ++pos){
01441                         m_miQuadRowIndexes[ i++] = pos->first;  
01442                 }
01443                 foundIdx.clear();       
01444                 return m_miQuadRowIndexes;
01445         }
01446         catch(const ErrorClass& eclass){
01447                 throw ErrorClass( eclass.errormsg);
01448         } 
01449 }//getQuadraticRowIndexes
01450 
01451 
01452 int OSInstance::getNumberOfQuadraticRowIndexes() {
01453         if(m_bQuadraticRowIndexesProcessed == false) getQuadraticRowIndexes();
01454         return m_iNumberOfQuadraticRowIndexes;
01455 }//getNumberOfQuadraticRowIndexes
01456 
01457 int* OSInstance::getNonlinearExpressionTreeIndexes(){
01458         if(m_bNonlinearExpressionTreeIndexesProcessed == true) return m_miNonlinearExpressionTreeIndexes;
01459         m_bNonlinearExpressionTreeIndexesProcessed = true;
01460         std::map<int, OSExpressionTree*> expTrees;
01461         expTrees = getAllNonlinearExpressionTrees();    
01462         std::map<int, OSExpressionTree*>::iterator pos;
01463         try{
01464                 // now put the term into an array
01465                 m_iNumberOfNonlinearExpressionTreeIndexes = expTrees.size();
01466                 m_miNonlinearExpressionTreeIndexes = new int[ m_iNumberOfNonlinearExpressionTreeIndexes ]       ;
01467                 int i = 0;
01468                 for(pos = expTrees.begin(); pos != expTrees.end(); ++pos){
01469                         m_miNonlinearExpressionTreeIndexes[ i++] = pos->first;  
01470                 }
01471                 expTrees.clear();       
01472                 return m_miNonlinearExpressionTreeIndexes;
01473         }
01474         catch(const ErrorClass& eclass){
01475                 throw ErrorClass( eclass.errormsg);
01476         } 
01477 }//getNonlinearExpressionTreeIndexes
01478 
01479 int OSInstance::getNumberOfNonlinearExpressionTreeIndexes() {
01480         if(m_bNonlinearExpressionTreeIndexesProcessed == false) getNonlinearExpressionTreeIndexes();
01481         return m_iNumberOfNonlinearExpressionTreeIndexes;
01482 }//getNumberOfNonlinearExpressionTreeIndexes
01483 
01484 
01485 
01486 int* OSInstance::getNonlinearExpressionTreeModIndexes(){
01487         if(m_bNonlinearExpressionTreeModIndexesProcessed == true) return m_miNonlinearExpressionTreeModIndexes;
01488         m_bNonlinearExpressionTreeModIndexesProcessed = true;
01489         std::map<int, OSExpressionTree*> expTrees;
01490         expTrees = getAllNonlinearExpressionTreesMod(); 
01491         std::map<int, OSExpressionTree*>::iterator pos;
01492         try{
01493                 // now put the term into an array
01494                 m_iNumberOfNonlinearExpressionTreeModIndexes = expTrees.size();
01495                 m_miNonlinearExpressionTreeModIndexes = new int[ m_iNumberOfNonlinearExpressionTreeModIndexes ] ;
01496                 int i = 0;
01497                 for(pos = expTrees.begin(); pos != expTrees.end(); ++pos){
01498                         m_miNonlinearExpressionTreeModIndexes[ i++] = pos->first;       
01499                 }
01500                 expTrees.clear();       
01501                 return m_miNonlinearExpressionTreeModIndexes;
01502         }
01503         catch(const ErrorClass& eclass){
01504                 throw ErrorClass( eclass.errormsg);
01505         } 
01506 }//getNonlinearExpressionTreeModIndexes
01507 
01508 int OSInstance::getNumberOfNonlinearExpressionTreeModIndexes() {
01509         if(m_bNonlinearExpressionTreeModIndexesProcessed == false) getNonlinearExpressionTreeModIndexes();
01510         return m_iNumberOfNonlinearExpressionTreeModIndexes;
01511 }//getNumberOfNonlinearExpressionTreeModIndexes
01512 
01513 
01514 int OSInstance::getNumberOfNonlinearConstraints(){
01515         if( m_bProcessExpressionTrees == false ){
01516                 getAllNonlinearExpressionTrees();
01517                 return m_iConstraintNumberNonlinear;
01518         }
01519         else return m_iConstraintNumberNonlinear;
01520 }//getNumberOfNonlinearConstraints
01521 
01522 int OSInstance::getNumberOfNonlinearObjectives(){
01523         if( m_bProcessExpressionTrees == false ){
01524                 getAllNonlinearExpressionTrees();
01525                 return m_iObjectiveNumberNonlinear;
01526         }
01527         else return m_iObjectiveNumberNonlinear;
01528 }//getNumberOfNonlinearObjectivess
01529 
01530 
01531 OSExpressionTree* OSInstance::getNonlinearExpressionTree(int rowIdx){
01532         if( m_bProcessExpressionTrees == false ){
01533                 getAllNonlinearExpressionTrees();
01534         } 
01535         if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()) return m_mapExpressionTrees[ rowIdx];
01536         else return NULL ;
01537         // check to make sure rowIdx has a nonlinear term and is in the map
01539         //std::map<int, OSExpressionTree*>::iterator pos;
01541         //      if(pos->first == rowIdx)return m_mapExpressionTrees[ rowIdx];
01542         //}
01543         // if we are rowIdx has no nonlinar terms so return a null
01544         //return NULL;
01545 }// getNonlinearExpressionTree for a specific index   
01546 
01547 
01548 OSExpressionTree* OSInstance::getNonlinearExpressionTreeMod(int rowIdx){
01549         if( m_bProcessExpressionTreesMod == false ){
01550                 getAllNonlinearExpressionTreesMod();
01551         } 
01552         if( m_mapExpressionTreesMod.find( rowIdx) != m_mapExpressionTreesMod.end()) return m_mapExpressionTreesMod[ rowIdx];
01553         else return NULL ;
01554         // check to make sure rowIdx has a nonlinear term and is in the map
01556         //std::map<int, OSExpressionTree*>::iterator pos;
01558         //      if(pos->first == rowIdx)return m_mapExpressionTrees[ rowIdx];
01559         //}
01560         // if we are rowIdx has no nonlinar terms so return a null
01561         //return NULL;
01562 }// getNonlinearExpressionTreeMod for a specific index 
01563 
01564 
01565 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeInPostfix( int rowIdx){
01566         //if( m_binitForAlgDiff == false) this->initForAlgDiff();
01567         if( m_bProcessExpressionTrees == false ) getAllNonlinearExpressionTrees();
01568         std::vector<OSnLNode*> postfixVec;
01569         try{
01570                 if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()){
01571                         OSExpressionTree* expTree = getNonlinearExpressionTree( rowIdx);
01572                         postfixVec = expTree->m_treeRoot->getPostfixFromExpressionTree();
01573                         
01574                 }  
01575                 else{
01576                         throw ErrorClass("Error in getNonlinearExpressionTreeInPostfix, rowIdx not valid");
01577                 }
01578                 return postfixVec;      
01579         }
01580         catch(const ErrorClass& eclass){
01581                 throw ErrorClass( eclass.errormsg);
01582         } 
01583 }//getNonlinearExpressionTreeInPostfix
01584 
01585 
01586 
01587 
01588 std::string OSInstance::getNonlinearExpressionTreeInInfix( int rowIdx_){
01589         if( m_binitForAlgDiff == false) this->initForAlgDiff();
01590         if( m_bProcessExpressionTrees == false ) getAllNonlinearExpressionTrees();
01591         std::string resultString;
01592         resultString = "";
01593         unsigned int i;
01594         unsigned int j;
01595         unsigned int n;
01596         ostringstream outStr;
01597         std::vector<OSnLNode*> postfixVec;
01598         int rowIdx = rowIdx_;
01599         OSnLNode *nlnode = NULL ;
01600         OSnLNodeNumber *nlnodeNum = NULL;
01601         OSnLNodeVariable *nlnodeVar = NULL;
01602         OSnLNodeSum *nlnodeSum = NULL;
01603         OSnLNodeProduct *nlnodeProduct = NULL;
01604         OSnLNodeMin *nlnodeMin = NULL;
01605         OSnLNodeMax *nlnodeMax = NULL;
01606         std::string tmp1 = "";
01607         std::string tmp2 = "";
01608         std::string tmp3 = "";
01609         std::stack<OSnLNode*> opStack;
01610         std::stack<std::string> tmpStack;
01611         std::stack<std::string> sumStack;
01612         std::stack<std::string> productStack;
01613         std::stack<std::string> minStack;
01614         std::stack<std::string> maxStack;
01615         
01616         
01617         try{
01618                 if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()){
01619                         //get the nodes and separate into operators and operands, for now
01620                         //only the number and variable nodes are operator nodes
01621                         
01622                         OSExpressionTree* exptree = this->getNonlinearExpressionTree( rowIdx);
01623                         if(exptree != NULL) {
01624                                 
01625                                 postfixVec = this->getNonlinearExpressionTreeInPostfix( rowIdx);
01626                                 n  = postfixVec.size();
01627                                 //put vector in reverse order
01628                                 for (i = 0 ; i < n; i++){
01629                                         nlnode =  postfixVec[ n - 1 - i];
01630                                         opStack.push( nlnode);
01631                                         //std::cout << postfixVec[ i]->snodeName << std::endl;
01632                                 }
01633                                 
01634                                 #if 0
01635                                 std::cout << std::endl << std::endl << std::endl << std::endl;
01636                                 n = opStack.size();
01637                                 for(i = 0; i < n; i++){
01638                                         std::cout  << opStack.top()->snodeName  << std::endl;
01639                                 }
01640                                 
01641                                 #endif
01642                                 n = opStack.size();
01643                                 for(i = 0; i < n; i++){
01644                                         //std::cout << "NUMBER OF NODES LEFT =  " << operatorVec.size() << std::endl;
01645                                         nlnode = opStack.top();
01646                                         //std::cout << "EVALUATING NODE: " << nlnode->snodeName << std::endl;
01647                                         switch (nlnode->inodeInt) {
01648                                                 case OS_NUMBER:
01649                                                         nlnodeNum = (OSnLNodeNumber*)nlnode;
01650                                                         tmpStack.push( os_dtoa_format(nlnodeNum->value) );
01651                                                         break;
01652                                                         
01653                                                 case OS_PI:
01654                                                         tmpStack.push( "PI" );
01655                                                         break;
01656                                                         
01657                                                 case OS_E:
01658                                                         tmpStack.push( "E" );
01659                                                         break;
01660                                                         
01661                                                 case OS_VARIABLE:
01662                                                         outStr.str("");
01663                                                         // handle a variable
01664                                                         nlnodeVar = (OSnLNodeVariable*)nlnode;
01665                                                         // see if the coefficient is specified
01666                                                         if( (nlnodeVar->coef > 1.0) ||  (nlnodeVar->coef < 1.0) ){
01667                                                                 outStr << "(";
01668                                                                 outStr <<  os_dtoa_format(nlnodeVar->coef);
01669                                                                 outStr << "*x_";
01670                                                                 outStr << nlnodeVar->idx;
01671                                                                 outStr << ")";
01672                                                                 tmpStack.push(outStr.str() );
01673                                                                 //std::cout << "WE JUST PUSHED " << outStr.str() << std::endl;
01674                                                         }else{
01675                                                                 outStr << "x_";
01676                                                                 outStr << nlnodeVar->idx;
01677                                                                 tmpStack.push(outStr.str() );
01678                                                                 //std::cout << "WE JUST PUSHED " << outStr.str() << std::endl;
01679                                                         }
01680                                                         break;
01681                                                         
01682                                                 case OS_PLUS :
01683                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing plus operater");
01684                                                         tmp1 = tmpStack.top();
01685                                                         tmpStack.pop();
01686                                                         tmp2 = tmpStack.top();
01687                                                         tmpStack.pop();
01688                                                         tmpStack.push("(" + tmp2 +  " + "  + tmp1 + ")");
01689                                                         break;
01690                                                         
01691                                                 case OS_SUM :
01692                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing sum operater");
01693                                                         //std::cout << "INSIDE SUM NODE " << std::endl;
01694                                                         nlnodeSum = (OSnLNodeSum*)nlnode;
01695                                                         outStr.str("");
01696                                                         for(j = 0; j < nlnodeSum->inumberOfChildren; j++){
01697                                                                 sumStack.push( tmpStack.top() );
01698                                                                 tmpStack.pop();
01699                                                         }
01700                                                         outStr << "(";
01701                                                         for(j = 0; j < nlnodeSum->inumberOfChildren; j++){
01702                                                                 outStr << sumStack.top();
01703                                                                 if (j < nlnodeSum->inumberOfChildren - 1) outStr << " + ";
01704                                                                 sumStack.pop();
01705                                                         }
01706                                                         outStr << ")";
01707                                                         tmpStack.push( outStr.str() );
01708                                                         //std::cout << outStr.str() << std::endl;
01709                                                         break;
01710                                                         
01711                                                 case OS_MINUS :
01712                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing minus operater");
01713                                                         tmp1 = tmpStack.top();
01714                                                         tmpStack.pop();
01715                                                         tmp2 = tmpStack.top();
01716                                                         tmpStack.pop();
01717                                                         tmpStack.push("(" + tmp2 +  " - "  + tmp1 + ")");
01718                                                         break;                                          
01719                                                 
01720                                                 case OS_NEGATE :
01721                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- -- Problem writing negate operater");
01722                                                         tmp1 = tmpStack.top();
01723                                                         tmpStack.pop();
01724                                                         tmpStack.push( "-"+ tmp1 );
01725                                                         break;
01726                                                         
01727                                                 case OS_TIMES :
01728                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing times operater");
01729                                                         tmp1 = tmpStack.top();
01730                                                         tmpStack.pop();
01731                                                         tmp2 = tmpStack.top();
01732                                                         tmpStack.pop();
01733                                                         tmpStack.push("(" + tmp2 +  "*"  + tmp1 + ")");
01734                                                         break;
01735                                                 
01736                                                 case OS_DIVIDE :
01737                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing divide operater");
01738                                                         tmp1 = tmpStack.top();
01739                                                         tmpStack.pop();
01740                                                         tmp2 = tmpStack.top();
01741                                                         tmpStack.pop();
01742                                                         tmpStack.push("(" + tmp2 +  " / "  + tmp1 + ")");
01743                                                         break;
01744 
01745                                                 case OS_POWER :
01746                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing power operater");
01747                                                         tmp1 = tmpStack.top();
01748                                                         tmpStack.pop();
01749                                                         tmp2 = tmpStack.top();
01750                                                         tmpStack.pop();
01751                                                         tmpStack.push("(" + tmp2 +  " ^ "  + tmp1 + ")");
01752                                                         break;
01753                                                         
01754                                                 case OS_ABS :
01755                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing abs operater");
01756                                                         tmp1 = tmpStack.top();
01757                                                         tmpStack.pop();
01758                                                         tmpStack.push( "abs( "+ tmp1  + ")");
01759                                                         break;
01760                                                         
01761                                                 case OS_ERF :
01762                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing erf operater");
01763                                                         tmp1 = tmpStack.top();
01764                                                         tmpStack.pop();
01765                                                         tmpStack.push( "erf( "+ tmp1  + ")");
01766                                                         break;
01767 
01768                                                 
01769                                                 case OS_SQUARE :
01770                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing square operater ");
01771                                                         tmp1 = tmpStack.top();
01772                                                         tmpStack.pop();
01773                                                         tmpStack.push( "("+ tmp1  + ")^2");
01774                                                         break;
01775                                                                                                         
01776                                                 case OS_LN :
01777                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing ln operater");
01778                                                         tmp1 = tmpStack.top();
01779                                                         tmpStack.pop();
01780                                                         tmpStack.push( "ln( "+ tmp1  + ")");
01781                                                         break;
01782                                                         
01783                                                 case OS_EXP :
01784                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing exp operater");
01785                                                         tmp1 = tmpStack.top();
01786                                                         tmpStack.pop();
01787                                                         tmpStack.push( "exp( "+ tmp1  + ")");
01788                                                         break;
01789                                                         
01790                                                 case OS_SIN :
01791                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing sin operater");
01792                                                         tmp1 = tmpStack.top();
01793                                                         tmpStack.pop();
01794                                                         tmpStack.push( "sin( "+ tmp1  + ")");
01795                                                         break;
01796                                                         
01797                                                 case OS_COS :
01798                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing cos operater ");
01799                                                         tmp1 = tmpStack.top();
01800                                                         tmpStack.pop();
01801                                                         tmpStack.push( "cos( "+ tmp1  + ")");
01802                                                         break;
01803                                                         
01804                                                 case OS_SQRT :
01805                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing cos operater ");
01806                                                         tmp1 = tmpStack.top();
01807                                                         tmpStack.pop();
01808                                                         tmpStack.push( "sqrt( "+ tmp1  + ")");
01809                                                         break;
01810                                                         
01811                                                 case OS_MIN :
01812                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing min operater");
01813                                                         //std::cout << "INSIDE Min NODE " << std::endl;
01814                                                         nlnodeMin = (OSnLNodeMin*)nlnode;
01815                                                         outStr.str("");
01816                                                         for(j = 0; j < nlnodeMin->inumberOfChildren; j++){
01817                                                                 minStack.push( tmpStack.top() );
01818                                                                 tmpStack.pop();
01819                                                         }
01820                                                         outStr << "min(";
01821                                                         for(j = 0; j < nlnodeMin->inumberOfChildren; j++){
01822                                                                 outStr << minStack.top();
01823                                                                 if (j < nlnodeMin->inumberOfChildren - 1) outStr << " , ";
01824                                                                 minStack.pop();
01825                                                         }
01826                                                         outStr << ")";
01827                                                         tmpStack.push( outStr.str() );
01828                                                         break;
01829                                                         
01830                                                 case OS_MAX :
01831                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing max operater");
01832                                                         //std::cout << "INSIDE Max NODE " << std::endl;
01833                                                         nlnodeMax = (OSnLNodeMax*)nlnode;
01834                                                         outStr.str("");
01835                                                         for(j = 0; j < nlnodeMax->inumberOfChildren; j++){
01836                                                                 maxStack.push( tmpStack.top() );
01837                                                                 tmpStack.pop();
01838                                                         }
01839                                                         outStr << "max(";
01840                                                         for(j = 0; j < nlnodeMax->inumberOfChildren; j++){
01841                                                                 outStr << maxStack.top();
01842                                                                 if (j < nlnodeMax->inumberOfChildren - 1) outStr << " , ";
01843                                                                 maxStack.pop();
01844                                                         }
01845                                                         outStr << ")";
01846                                                         tmpStack.push( outStr.str() );
01847                                                         break;
01848                                                         
01849                                                 case OS_IF :
01850                                                 
01851                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing if operater ");
01852                                                         if(nlnode->inumberOfChildren != 3)throw  ErrorClass("The if node must have three children");
01853                                                         tmp1 = tmpStack.top();
01854                                                         tmpStack.pop();
01855                                                         tmp2 = tmpStack.top();
01856                                                         tmpStack.pop();
01857                                                         tmp3 = tmpStack.top();
01858                                                         tmpStack.pop();
01859                                                         tmpStack.push( "if(" + tmp3 + "," + tmp2 + "," + tmp1 +")" );
01860                                                         break;
01861                                                         
01862                                                 case OS_PRODUCT :
01863                                                         if( tmpStack.size() < nlnode->inumberOfChildren) throw  ErrorClass("There is an error in the OSExpression Tree -- Problem writing product operater");
01864                                                         //std::cout << "INSIDE Product NODE " << std::endl;
01865                                                         nlnodeProduct = (OSnLNodeProduct*)nlnode;
01866                                                         outStr.str("");
01867                                                         for(j = 0; j < nlnodeProduct->inumberOfChildren; j++){
01868                                                                 productStack.push( tmpStack.top() );
01869                                                                 tmpStack.pop();
01870                                                         }
01871                                                         outStr << "(";
01872                                                         for(j = 0; j < nlnodeProduct->inumberOfChildren; j++){
01873                                                                 outStr << productStack.top();
01874                                                                 if (j < nlnodeProduct->inumberOfChildren - 1) outStr << " * ";
01875                                                                 productStack.pop();
01876                                                         }
01877                                                         outStr << ")";
01878                                                         tmpStack.push( outStr.str() );
01879                                                         //std::cout << outStr.str() << std::endl;
01880                                                         break;
01881                                                         
01882                                                 default:
01883                                                          throw  ErrorClass("operator " + nlnode->snodeName + " not supported");
01884                                                 break;
01885                                   }
01886                                         opStack.pop();
01887                                 }
01888                                 postfixVec.clear();
01889                                 if(tmpStack.size() != 1) throw ErrorClass( "There is an error in the OSExpression Tree -- stack size should be 1 at end");
01890                                 resultString = tmpStack.top();
01891                                 //std::cout << resultString << std::endl;
01892                                 tmpStack.pop();
01893                                 
01894                                 return resultString;
01895 
01896                         }
01897                         else{
01898                                 //throw ErrorClass("Error in getNonlinearExpressionTreeInInfix, there is no expression tree for this index");
01899                                 return "";
01900                         }
01901                         
01902                 }  
01903                 else{
01904                         throw ErrorClass("Error in getNonlinearExpressionTreeInInfix, rowIdx not valid");
01905                 }
01906                 return resultString;    
01907         }
01908         catch(const ErrorClass& eclass){
01909                 throw ErrorClass( eclass.errormsg);
01910         } 
01911 }//getNonlinearExpressionTreeInInfix
01912 
01913 std::string OSInstance::printModel( ){
01914         std::string resultString = "";
01915         ostringstream outStr;
01916         outStr << "";
01917         int numCon;
01918         int numObj;
01919         int numVar;
01920         int i;
01921         numCon = this->getConstraintNumber();
01922         numObj = this->getObjectiveNumber();
01923         numVar = this->getVariableNumber();
01924         this->initForAlgDiff( );
01925         outStr << std::endl;
01926         outStr << "Objectives:" << std::endl;           
01927         for(i = 0; i < numObj; i++){
01928                 if( this->getObjectiveMaxOrMins()[0] == "min"){
01929                         outStr << "min ";
01930                 }
01931                 else{
01932                         outStr << "max ";
01933                 }
01934                 outStr << this->printModel( i - numObj);
01935         }       
01936         outStr << std::endl;
01937         outStr << "Constraints:" << std::endl;          
01938         for(i = 0; i < numCon; i++){
01939                 outStr << this->printModel( i);
01940         }
01941         
01942         outStr << std::endl;
01943         outStr << "Variables:" << std::endl;
01944         if(m_bProcessVariables != true ) this->processVariables();
01945         for(i = 0; i < numVar; i++){
01946                 outStr << "x_";
01947                 outStr << i;
01948                 outStr << " Type = " ;
01949                 outStr <<  m_mcVariableTypes[i];
01950                 outStr << "  Lower Bound =  ";
01951                 outStr << os_dtoa_format( m_mdVariableLowerBounds[i])  ;
01952                 outStr << "  Upper Bound =  ";
01953                 outStr << os_dtoa_format( m_mdVariableUpperBounds[i])  ;
01954                 outStr << std::endl;
01955         }       
01956         
01957         // if model was originaly in column matrix form we need to delete the new
01958         // matrix stored by row
01959         
01960         //if( this->instanceData->linearConstraintCoefficients != NULL && this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
01961         //      if(m_bColumnMajor == true){// we now have a duplication
01962         //              delete m_linearConstraintCoefficientsInRowMajor;
01963         //      }
01964         //}
01965         return outStr.str() ;
01966 }//printModel( )
01967 
01968 
01969 std::string OSInstance::printModel(int rowIdx ){
01970         std::string resultString = "";
01971         ostringstream outStr;
01972         outStr << "";
01973         //loop over the consraints first;
01974         int j;
01975         int row_nonz = 0;
01976         int obj_nonz = 0;
01977         int varIdx = 0;
01978         bool addedLinearTerm = false;
01979         
01980         // initialize all of the necessary nonlinear stuff
01981         this->initForAlgDiff( );        
01982         
01983         if( rowIdx >= 0){
01984                 if( m_bProcessConstraints != true ) this->processConstraints() ;
01985                 if( m_mdConstraintLowerBounds[ rowIdx] >  -OSDBL_MAX){
01986                         if(m_mdConstraintLowerBounds[ rowIdx] < m_mdConstraintUpperBounds[ rowIdx]){
01987                                 outStr << os_dtoa_format( m_mdConstraintLowerBounds[ rowIdx] );
01988                                 outStr << " <= ";
01989                         }
01990                 } 
01991                 //
01992                 if(this->instanceData->linearConstraintCoefficients != NULL && this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
01993                         if(m_linearConstraintCoefficientsInRowMajor == NULL) 
01994                                 m_linearConstraintCoefficientsInRowMajor = this->getLinearConstraintCoefficientsInRowMajor();
01995                         row_nonz = m_linearConstraintCoefficientsInRowMajor->starts[ rowIdx + 1] - m_linearConstraintCoefficientsInRowMajor->starts[ rowIdx];
01996                         
01997                         for(j = 0; j < row_nonz; j++){
01998                                 varIdx =  m_linearConstraintCoefficientsInRowMajor->indexes[ m_linearConstraintCoefficientsInRowMajor->starts[ rowIdx]  + j];
01999                                 
02000                                 if(m_bSparseJacobianCalculated == false  ||  (m_mapExpressionTreesMod.find( rowIdx) == m_mapExpressionTreesMod.end() ) ||
02001                                         ( (*m_mapExpressionTreesMod[ rowIdx]->mapVarIdx).find( varIdx) == (*m_mapExpressionTreesMod[ rowIdx]->mapVarIdx).end()) ){                                      
02002                                         outStr << os_dtoa_format( m_linearConstraintCoefficientsInRowMajor->values[ m_linearConstraintCoefficientsInRowMajor->starts[ rowIdx]  + j] );
02003                                         outStr << "*";
02004                                         outStr << "x_";
02005                                         outStr << varIdx;
02006                                         if( j < row_nonz - 1) outStr << " + ";
02007                                         addedLinearTerm = true;
02008                                 }
02009                         }
02010                 }
02011         }else{// process an objective function
02012                 if(m_bProcessObjectives != true ) this->processObjectives() ;
02013                 int obj_idx =  -rowIdx - 1;
02014                 obj_nonz = m_miNumberOfObjCoef[ obj_idx];
02015                 for(j = 0; j < obj_nonz; j++){
02016                         outStr << os_dtoa_format( m_mObjectiveCoefficients[obj_idx]->values[j] );
02017                         outStr << "*";
02018                         outStr << "x_";
02019                         outStr << m_mObjectiveCoefficients[obj_idx]->indexes[j] ;
02020                         if( j < obj_nonz - 1) outStr << " + ";          
02021                 }
02022         }
02023         if( this->getNonlinearExpressionTree( rowIdx) != NULL){
02024                 if( (addedLinearTerm == true)  || (obj_nonz > 0) ) outStr << " + " ;
02025                 outStr << getNonlinearExpressionTreeInInfix( rowIdx);
02026                 //outStr << ")";
02027         }
02028         
02029         
02030         if( rowIdx >= 0){
02031                 if( m_bProcessConstraints != true ) this->processConstraints() ;
02032                 if( m_mdConstraintUpperBounds[ rowIdx] <  OSDBL_MAX){
02033                         if(m_mdConstraintLowerBounds[ rowIdx] < m_mdConstraintUpperBounds[ rowIdx]){
02034                                 outStr << " <= ";
02035                          }else{
02036                                 outStr << " = ";
02037                          }
02038                         outStr << os_dtoa_format( m_mdConstraintUpperBounds[ rowIdx] ); 
02039                         outStr << "  ";
02040                         outStr << this->getConstraintNames()[ rowIdx] ;
02041                 }
02042         }
02043         outStr << std::endl;    
02044         resultString = outStr.str();
02045         return resultString;
02046 }//printModel( rowIdx )
02047 
02048 
02049 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeModInPostfix( int rowIdx){
02050         //if( m_binitForAlgDiff == false) this->initForAlgDiff();
02051         if( m_bProcessExpressionTreesMod == false ) getAllNonlinearExpressionTreesMod();
02052         std::vector<OSnLNode*> postfixVec;
02053         try{
02054                 if( m_mapExpressionTreesMod.find( rowIdx) != m_mapExpressionTreesMod.end()){
02055                         OSExpressionTree* expTree = getNonlinearExpressionTreeMod( rowIdx);
02056                         postfixVec = expTree->m_treeRoot->getPostfixFromExpressionTree();
02057                         
02058                 }  
02059                 else{
02060                         throw ErrorClass("Error in getNonlinearExpressionTreeModInPostfix, rowIdx not valid");
02061                 }
02062                 return postfixVec;      
02063         }
02064         catch(const ErrorClass& eclass){
02065                 throw ErrorClass( eclass.errormsg);
02066         } 
02067 }//getNonlinearExpressionTreeModInPostfix
02068 
02069 
02070 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeInPrefix( int rowIdx){
02071         //if( m_binitForAlgDiff == false) this->initForAlgDiff();
02072         if( m_bProcessExpressionTrees == false ) getAllNonlinearExpressionTrees();
02073         std::vector<OSnLNode*> prefixVec;
02074         try{
02075                 if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()){
02076                         OSExpressionTree* expTree = getNonlinearExpressionTree( rowIdx);
02077                         prefixVec = expTree->m_treeRoot->getPrefixFromExpressionTree();
02078                         
02079                 }  
02080                 else{
02081                         throw ErrorClass("Error in getNonlinearExpressionTreeInPrefix, rowIdx not valid");
02082                 }
02083                 return prefixVec;       
02084         }
02085         catch(const ErrorClass& eclass){
02086                 throw ErrorClass( eclass.errormsg);
02087         } 
02088 }//getNonlinearExpressionTreeInPrefix
02089 
02090 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeModInPrefix( int rowIdx){
02091         //if( m_binitForAlgDiff == false) this->initForAlgDiff();
02092         if( m_bProcessExpressionTreesMod == false ) getAllNonlinearExpressionTreesMod();
02093         std::vector<OSnLNode*> prefixVec;
02094         try{
02095                 if( m_mapExpressionTreesMod.find( rowIdx) != m_mapExpressionTreesMod.end()){
02096                         OSExpressionTree* expTree = getNonlinearExpressionTreeMod( rowIdx);
02097                         prefixVec = expTree->m_treeRoot->getPrefixFromExpressionTree();
02098                         
02099                 }  
02100                 else{
02101                         throw ErrorClass("Error in getNonlinearExpressionTreeInPrefix, rowIdx not valid");
02102                 }
02103                 return prefixVec;       
02104         }
02105         catch(const ErrorClass& eclass){
02106                 throw ErrorClass( eclass.errormsg);
02107         } 
02108 }//getNonlinearExpressionTreeInPrefix
02109 
02110 std::map<int, OSExpressionTree*> OSInstance::getAllNonlinearExpressionTrees(){
02111         //if( m_binitForAlgDiff == false) this->initForAlgDiff();
02112         if(m_bProcessExpressionTrees == true) return m_mapExpressionTrees;
02113         std::map<int, int> foundIdx;
02114         std::map<int, int>::iterator pos;
02115         OSnLNodePlus *nlNodePlus;
02116         OSExpressionTree *expTree;  
02117         m_iObjectiveNumberNonlinear = 0;   
02118         m_iConstraintNumberNonlinear = 0;    
02119         int i;   
02120         // important -- tell the nl nodes not to destroy the OSExpression Objects
02121         if( instanceData->nonlinearExpressions->numberOfNonlinearExpressions > 0 && instanceData->nonlinearExpressions->nl != NULL){
02122                 for( i = 0; i < instanceData->nonlinearExpressions->numberOfNonlinearExpressions; i++){
02123                         instanceData->nonlinearExpressions->nl[i]->m_bDeleteExpressionTree = false;
02124                 }
02125         }
02126         int index;  
02127         // kipp -- what should we return if instanceData->nonlinearExpressions->numberOfNonlinearExpressions is zero
02128         for(i = 0; i < instanceData->nonlinearExpressions->numberOfNonlinearExpressions; i++){
02129                 index = instanceData->nonlinearExpressions->nl[ i]->idx;
02130                 if(foundIdx.find( index) != foundIdx.end() ){  
02131                 //if(foundIdx[ index] > 0 ){ 
02132                         //std::cout << "OLD INDEX FOUND " << index << std::endl;
02133                         //std::cout << "foundIdx[ index] " << index << std::endl;
02134                         // found an existing index
02135                         // important -- at this time m_mapExpressionTrees[ index] points to 
02136                         // the last OSExpressionTree with this index, it does not point to the 
02137                         // the just found OSExpressionTree with this index
02138                         nlNodePlus = new OSnLNodePlus();
02139                         //expTree = new OSExpressionTree(); 
02140                         expTree =  instanceData->nonlinearExpressions->nl[ i]->osExpressionTree;
02141                         // set left child to old index and right child to new one 
02142                         nlNodePlus->m_mChildren[ 0] = m_mapExpressionTrees[ index]->m_treeRoot;
02143                         nlNodePlus->m_mChildren[ 1] = instanceData->nonlinearExpressions->nl[ i]->osExpressionTree->m_treeRoot;
02144                         // we must delete the Expression tree corresponding to the old index value but not the nl nodes
02145                         instanceData->nonlinearExpressions->nl[ foundIdx[ index]  ]->m_bDeleteExpressionTree = true;
02146                         instanceData->nonlinearExpressions->nl[ foundIdx[ index]  ]->osExpressionTree->bDestroyNlNodes = false;
02147                         //point to the new expression tree
02148                         m_mapExpressionTrees[ index] = expTree;
02149                         m_mapExpressionTrees[ index]->m_treeRoot = nlNodePlus;
02150                         foundIdx[ index] = i;
02151                 }
02152                 else{  
02153                         // we have a new index
02154                         m_mapExpressionTrees[ index] = instanceData->nonlinearExpressions->nl[ i]->osExpressionTree;
02155                         m_mapExpressionTrees[ index]->m_treeRoot = instanceData->nonlinearExpressions->nl[ i]->osExpressionTree->m_treeRoot;
02156                         foundIdx[ index] = i;
02157                 }
02158                 //foundIdx[ index]++;    
02159         }
02160         // count the number of constraints and objective functions with nonlinear terms.
02161         for(pos = foundIdx.begin(); pos != foundIdx.end(); ++pos){
02162                 if(pos->first == -1) {
02163                         m_iObjectiveNumberNonlinear++; 
02164                 }
02165                 else m_iConstraintNumberNonlinear++;
02166         }
02167         foundIdx.clear();
02168         m_bProcessExpressionTrees = true;
02169         return m_mapExpressionTrees;
02170 }// getAllNonlinearExpressionTrees
02171 
02172 std::map<int, OSExpressionTree*> OSInstance::getAllNonlinearExpressionTreesMod(){
02173         if( m_bProcessExpressionTreesMod == true ) return m_mapExpressionTreesMod;
02174         m_bProcessExpressionTreesMod = true;
02175         // make sure we have the modified map available
02176         if( m_bNonLinearStructuresInitialized == false) initializeNonLinearStructures( );
02177         return m_mapExpressionTreesMod;
02178 }// getAllNonlinearExpressionTreesMod
02179 
02180 
02186 std::string OSInstance::getTimeDomainFormat()
02187 {       if (instanceData->timeDomain == NULL)
02188                 return "";
02189         if (instanceData->timeDomain->interval != NULL)
02190                 return "interval";
02191         if (instanceData->timeDomain->stages != NULL)
02192                 return "stages";
02193         return "";
02194 } 
02195 
02201 int OSInstance::getTimeDomainStageNumber()
02202 {       if (instanceData->timeDomain == NULL)
02203                 return 1;
02204         if (instanceData->timeDomain->interval != NULL)
02205                 ; //throw an error
02206         if (instanceData->timeDomain->stages == NULL)
02207                 return 1;
02208         return instanceData->timeDomain->stages->numberOfStages;
02209 } 
02210         
02216 std::string* OSInstance::getTimeDomainStageNames()
02217 {       if (instanceData->timeDomain == NULL)
02218                 return NULL;
02219         if (instanceData->timeDomain->interval != NULL)
02220                 return NULL; //throw an error
02221         if (instanceData->timeDomain->stages == NULL)
02222                 return NULL;
02223         if (m_msTimeDomainStageNames != NULL)
02224                 delete [] m_msTimeDomainStageNames;
02225         if (instanceData->timeDomain->stages->numberOfStages == 0)
02226                 return NULL;
02227         m_msTimeDomainStageNames = new std::string[instanceData->timeDomain->stages->numberOfStages]; 
02228         for (int i = 0; i < instanceData->timeDomain->stages->numberOfStages; i++)
02229                 m_msTimeDomainStageNames[i] = instanceData->timeDomain->stages->stage[i]->name;
02230         return m_msTimeDomainStageNames;
02231 } 
02232         
02238 int* OSInstance::getTimeDomainStageNumberOfVariables()
02239 {       if (instanceData->timeDomain == NULL)
02240                 return NULL;
02241         if (instanceData->timeDomain->interval != NULL)
02242                 return NULL; //throw an error
02243         if (instanceData->timeDomain->stages == NULL)
02244                 return NULL;
02245         if (m_miTimeDomainStageVariableNumber != NULL) 
02246                 delete [] m_miTimeDomainStageVariableNumber;
02247         if (instanceData->timeDomain->stages->numberOfStages == 0)
02248                 return NULL;
02249         m_miTimeDomainStageVariableNumber = new int[instanceData->timeDomain->stages->numberOfStages];
02250         for (int i = 0; i < instanceData->timeDomain->stages->numberOfStages; i++) 
02251                 m_miTimeDomainStageVariableNumber[i] = instanceData->timeDomain->stages->stage[i]->variables->numberOfVariables;
02252         return m_miTimeDomainStageVariableNumber;
02253 } 
02254         
02260 int* OSInstance::getTimeDomainStageNumberOfConstraints()
02261 {       if (instanceData->timeDomain == NULL)
02262                 return NULL;
02263         if (instanceData->timeDomain->interval != NULL)
02264                 return NULL; //throw an error
02265         if (instanceData->timeDomain->stages == NULL)
02266                 return NULL;
02267         if (m_miTimeDomainStageConstraintNumber != NULL){
02268         }
02269         
02270         
02271                 delete [] m_miTimeDomainStageConstraintNumber;
02272         if (instanceData->timeDomain->stages->numberOfStages == 0)
02273                 return NULL;
02274         m_miTimeDomainStageConstraintNumber = new int[instanceData->timeDomain->stages->numberOfStages];
02275         for (int i = 0; i < instanceData->timeDomain->stages->numberOfStages; i++)
02276                 m_miTimeDomainStageConstraintNumber[i] = instanceData->timeDomain->stages->stage[i]->constraints->numberOfConstraints;
02277         return m_miTimeDomainStageConstraintNumber;
02278 } 
02279 
02285 int* OSInstance::getTimeDomainStageNumberOfObjectives()
02286 {       if (instanceData->timeDomain == NULL)
02287                 return NULL;
02288         if (instanceData->timeDomain->interval != NULL)
02289                 return NULL; //throw an error
02290         if (instanceData->timeDomain->stages == NULL)
02291                 return NULL;
02292         if (m_miTimeDomainStageObjectiveNumber != NULL)
02293                 delete [] m_miTimeDomainStageObjectiveNumber;
02294         if (instanceData->timeDomain->stages->numberOfStages == 0)
02295                 return NULL;
02296         m_miTimeDomainStageObjectiveNumber = new int[instanceData->timeDomain->stages->numberOfStages];
02297         for (int i = 0; i < instanceData->timeDomain->stages->numberOfStages; i++)
02298                 m_miTimeDomainStageObjectiveNumber[i] = instanceData->timeDomain->stages->stage[i]->objectives->numberOfObjectives;
02299         return m_miTimeDomainStageObjectiveNumber;
02300 } 
02301         
02307 int** OSInstance::getTimeDomainStageVarList()
02308 {       if (instanceData->timeDomain == NULL)
02309                 return NULL;
02310         if (instanceData->timeDomain->interval != NULL)
02311                 return NULL; //throw an error
02312         if (instanceData->timeDomain->stages == NULL)
02313                 return NULL;
02314         if (m_miTimeDomainStageVariableNumber == NULL)
02315                 return NULL;
02316         if (m_mmiTimeDomainStageVarList != NULL){
02317                 
02318                 
02319                 for (int i = 0; i < m_iNumberOfTimeStages; i ++) 
02320                         delete[] m_mmiTimeDomainStageVarList[i];
02321                 delete[] m_mmiTimeDomainStageVarList;
02322                 m_mmiTimeDomainStageVarList = NULL;
02323         }
02324         
02325         
02326                 //delete [] m_mmiTimeDomainStageVarList;
02327         if (instanceData->timeDomain->stages->numberOfStages == 0)
02328                 return NULL;
02329         m_iNumberOfTimeStages = instanceData->timeDomain->stages->numberOfStages;
02330         m_mmiTimeDomainStageVarList = new int*[instanceData->timeDomain->stages->numberOfStages];
02331         int timeDomainStageNumberVar;
02332         for (int i = 0; i < instanceData->timeDomain->stages->numberOfStages; i++)
02333         {       
02334                 timeDomainStageNumberVar = m_miTimeDomainStageVariableNumber[i];
02335                 m_mmiTimeDomainStageVarList[i] = new int[ timeDomainStageNumberVar ];
02336                 if (instanceData->timeDomain->stages->stage[i]->variables->startIdx == -1)
02337                         for (int j = 0; j < m_miTimeDomainStageVariableNumber[i]; j++)
02338                                 m_mmiTimeDomainStageVarList[i][j] = instanceData->timeDomain->stages->stage[i]->variables->var[j]->idx;
02339                 else
02340                         for (int j = 0; j < m_miTimeDomainStageVariableNumber[i]; j++)
02341                                 m_mmiTimeDomainStageVarList[i][j] = instanceData->timeDomain->stages->stage[i]->variables->startIdx + j;
02342         }
02343         return m_mmiTimeDomainStageVarList;
02344 } 
02345         
02351 int** OSInstance::getTimeDomainStageConList()
02352 {       if (instanceData->timeDomain == NULL)
02353                 return NULL;
02354         if (instanceData->timeDomain->interval != NULL)
02355                 return NULL; //throw an error
02356         if (instanceData->timeDomain->stages == NULL)
02357                 return NULL;
02358         if (m_miTimeDomainStageConstraintNumber == NULL)
02359                 return NULL;
02360         if (m_mmiTimeDomainStageConList != NULL){
02361                 
02362                 for (int i = 0; i < m_iNumberOfTimeStages; i ++) 
02363                         delete[] m_mmiTimeDomainStageConList[i];
02364                 delete[] m_mmiTimeDomainStageConList;
02365                 m_mmiTimeDomainStageConList = NULL;             
02366                 
02367         }
02368         
02369         
02370                 
02371                 
02372                 
02373                 
02374                 
02375         if (instanceData->timeDomain->stages->numberOfStages == 0)
02376                 return NULL;
02377         m_iNumberOfTimeStages = instanceData->timeDomain->stages->numberOfStages;
02378         m_mmiTimeDomainStageConList = new int*[instanceData->timeDomain->stages->numberOfStages];
02379         int numTimeDomainStageCon;
02380         for (int i = 0; i < instanceData->timeDomain->stages->numberOfStages; i++)
02381         {       
02382                 numTimeDomainStageCon = m_miTimeDomainStageConstraintNumber[i];
02383                 m_mmiTimeDomainStageConList[i] = new int[ numTimeDomainStageCon];
02384                 if (instanceData->timeDomain->stages->stage[i]->constraints->startIdx == -1)
02385                         for (int j = 0; j < m_miTimeDomainStageConstraintNumber[i]; j++)
02386                                 m_mmiTimeDomainStageConList[i][j] = instanceData->timeDomain->stages->stage[i]->constraints->con[j]->idx;
02387                 else
02388                         for (int j = 0; j < m_miTimeDomainStageConstraintNumber[i]; j++)
02389                                 m_mmiTimeDomainStageConList[i][j] = instanceData->timeDomain->stages->stage[i]->constraints->startIdx + j;
02390         }
02391         return m_mmiTimeDomainStageConList;
02392 } 
02398 int** OSInstance::getTimeDomainStageObjList()
02399 {       if (instanceData->timeDomain == NULL)
02400                 return NULL;
02401         if (instanceData->timeDomain->interval != NULL)
02402                 return NULL; //throw an error
02403         if (instanceData->timeDomain->stages == NULL)
02404                 return NULL;
02405         if (m_miTimeDomainStageObjectiveNumber == NULL)
02406                 return NULL;
02407         if (m_mmiTimeDomainStageObjList != NULL){
02408 
02409                 for (int i = 0; i < m_iNumberOfTimeStages; i ++) 
02410                         delete[] m_mmiTimeDomainStageObjList[i];
02411                 delete[] m_mmiTimeDomainStageObjList;
02412                 m_mmiTimeDomainStageObjList = NULL;
02413                 
02414         }
02415         if (instanceData->timeDomain->stages->numberOfStages == 0)
02416                 return NULL;
02417         m_iNumberOfTimeStages = instanceData->timeDomain->stages->numberOfStages;
02418         m_mmiTimeDomainStageObjList = new int*[instanceData->timeDomain->stages->numberOfStages];
02419         int numTimeDomainStageObjNum;
02420         for (int i = 0; i < instanceData->timeDomain->stages->numberOfStages; i++)
02421         {       
02422                 numTimeDomainStageObjNum = m_miTimeDomainStageObjectiveNumber[i];
02423                 m_mmiTimeDomainStageObjList[i] = new int[ numTimeDomainStageObjNum];
02424                 if (instanceData->timeDomain->stages->stage[i]->objectives->startIdx == 0)
02425                         for (int j = 0; j < m_miTimeDomainStageObjectiveNumber[i]; j++)
02426                                 m_mmiTimeDomainStageObjList[i][j] = instanceData->timeDomain->stages->stage[i]->objectives->obj[j]->idx;
02427                 else
02428                         for (int j = 0; j < m_miTimeDomainStageObjectiveNumber[i]; j++)
02429                                 m_mmiTimeDomainStageObjList[i][j] = instanceData->timeDomain->stages->stage[i]->objectives->startIdx - j;
02430         }
02431         return m_mmiTimeDomainStageObjList;
02432 } 
02433 
02439 double OSInstance::getTimeDomainIntervalStart()
02440 {       if (instanceData->timeDomain == NULL)
02441                 return 0.0;
02442         if (instanceData->timeDomain->stages != NULL)
02443                 return 0.0; //throw an error
02444         if (instanceData->timeDomain->interval == NULL)
02445                 return 0.0;
02446         return instanceData->timeDomain->interval->start;
02447 } 
02448 
02454 double OSInstance::getTimeDomainIntervalHorizon()
02455 {       if (instanceData->timeDomain == NULL)
02456                 return 0.0;
02457         if (instanceData->timeDomain->stages != NULL)
02458                 return 0.0; //throw an error
02459         if (instanceData->timeDomain->interval == NULL)
02460                 return 0.0;
02461         return instanceData->timeDomain->interval->horizon;
02462 } 
02463 
02464 
02465 
02466 
02467 // the set() methods
02468 
02469 bool OSInstance::setInstanceSource(string source){
02470         instanceHeader->source = source;
02471         return true;
02472 }//setInstanceSource
02473 
02474 bool OSInstance::setInstanceDescription(string description){
02475         instanceHeader->description = description;
02476         return true;
02477 }//setInstanceDescription
02478 
02479 
02480 bool OSInstance::setInstanceName(string name){
02481          instanceHeader->name = name;
02482          return true;
02483 }//setInstanceName
02484 
02485 
02486 bool OSInstance::setVariableNumber(int number){
02487         // this method assume osinstance->instanceData->variables is not null
02488         if(number <= 0) return false;
02489         //if(instanceData->variables->numberOfVariables != -1  && instanceData->variables->numberOfVariables != number){
02490         //      delete[] instanceData->variables->var;
02491         //      instanceData->variables->var = NULL;
02492         //} 
02493         if(instanceData->variables == NULL) instanceData->variables = new Variables();
02494         instanceData->variables->numberOfVariables = number;
02495         if(instanceData->variables->var == NULL){
02496                 instanceData->variables->var = new Variable*[number];
02497         }
02498         return true;
02499 }//setVariableNumber
02500 
02501 
02502 bool OSInstance::addVariable(int index, string name, double lowerBound, double upperBound, char type){
02503         instanceData->variables->var[index] = new Variable();
02504         if(index < 0 || instanceData->variables->numberOfVariables <= 0 || index >= instanceData->variables->numberOfVariables) return false;
02505         instanceData->variables->var[index]->name = name;
02506         instanceData->variables->var[index]->lb = lowerBound;
02507         instanceData->variables->var[index]->ub = upperBound;
02508         instanceData->variables->var[index]->type = type;
02509         //if(init != OSNAN) instanceData->variables->var[index]->init = init;
02510         //instanceData->variables->var[index]->initString = initString;
02511         return true;
02512 }//addVariable
02513 
02514 
02515 bool OSInstance::setVariables(int number, string *names, double *lowerBounds, 
02516         double *upperBounds, char *types){
02517         if(number <= 0) return false;
02518         try{
02519                 if(instanceData->variables == NULL){
02520                         throw ErrorClass("There is no variables object");
02521                 }
02522                 if(instanceData->variables->numberOfVariables != number){
02523                         throw ErrorClass("input number of variables not equal to number in class");
02524                 }
02525                 //instanceData->variables->var = new Variable*[number];
02526                 int i;
02527                 for(i = 0; i < number; i++){
02528                         instanceData->variables->var[ i] = new Variable();
02529                 }
02530                 if(names  != NULL){
02531                         for(i = 0; i < number; i++) instanceData->variables->var[i]->name = names[i];                           
02532                 }
02533                 if(lowerBounds != NULL){
02534                         for(i = 0; i < number; i++){
02535                                 instanceData->variables->var[i]->lb = lowerBounds[i];  
02536                         }
02537                 }
02538                 if(upperBounds != NULL){
02539                         for(i = 0; i < number; i++){
02540                                 instanceData->variables->var[i]->ub = upperBounds[i]; 
02541                         }
02542                 }
02543                 if(types != NULL){
02544                         for(i = 0; i < number; i++){
02545                                 instanceData->variables->var[i]->type = types[i];
02546                                 if(types[i] != 'C' && types[i] != 'B' && types[i] != 'I' && types[i] != 'S') types[i] = 'C';
02547                         } 
02548                 }
02549                 /*
02550                 if(inits != NULL){
02551                         for(i = 0; i < number; i++) instanceData->variables->var[i]->init = inits[i];                           
02552                 }
02553                 if(initsString != NULL){
02554                         for(i = 0; i < number; i++) instanceData->variables->var[i]->initString = initsString[i];                       
02555                 }
02556                 */
02557                 return true;
02558         }
02559         catch(const ErrorClass& eclass){
02560                 throw ErrorClass(  eclass.errormsg); 
02561         }
02562 }//setVariables
02563 
02564 // begin checking again with Jun Ma
02565 
02566 bool OSInstance::setObjectiveNumber(int number){
02567         if(number < 0) return false;
02568         if(instanceData->objectives == NULL) instanceData->objectives = new Objectives();
02569         if(number == 0){
02570                 instanceData->objectives->numberOfObjectives = 0;
02571                 instanceData->objectives->obj = 0;
02572                 return true;
02573         }
02574         instanceData->objectives->numberOfObjectives = number;
02575         instanceData->objectives->obj = new Objective*[number];                         
02576         return true;
02577 }//setObjectiveNumber
02578 
02579 bool OSInstance::addObjective(int index, string name, string maxOrMin, double constant, double weight, SparseVector *objectiveCoefficients){
02580         if(index >= 0 || instanceData->objectives->numberOfObjectives <= 0 || abs(index) > instanceData->objectives->numberOfObjectives) return false;
02581         int arrayIndex = abs(index) -1;
02582         if(instanceData->objectives->obj == NULL) return false;
02583         instanceData->objectives->obj[arrayIndex] = new Objective();
02584         instanceData->objectives->obj[arrayIndex]->name = name;
02585         if( (maxOrMin != "max") && (maxOrMin != "min") ) return false;
02586         else instanceData->objectives->obj[arrayIndex]->maxOrMin = maxOrMin;
02587         instanceData->objectives->obj[arrayIndex]->constant = constant;
02588         instanceData->objectives->obj[arrayIndex]->weight = weight;
02589         int n = objectiveCoefficients->number;                  
02590         instanceData->objectives->obj[arrayIndex]->numberOfObjCoef = n;
02591         if(n == 0){
02592                 instanceData->objectives->obj[arrayIndex]->coef = 0;
02593         }
02594         else{
02595                 int i = 0;
02596                 instanceData->objectives->obj[arrayIndex]->coef = new ObjCoef*[n];
02597                 for(i = 0; i < n; i++) instanceData->objectives->obj[arrayIndex]->coef[i] = new ObjCoef();
02598                 for(i = 0; i < n; i++){
02599                         instanceData->objectives->obj[arrayIndex]->coef[i]->idx = objectiveCoefficients->indexes[i];
02600                         instanceData->objectives->obj[arrayIndex]->coef[i]->value = objectiveCoefficients->values[i];   
02601                 }                       
02602         }
02603         return true;
02604 }//addObjective
02605 
02606 bool OSInstance::setObjectives(int number, string *names, string *maxOrMins, double *constants, double *weights, SparseVector **objectiveCoefficients){
02607         if(number < 0) return false;
02608         try{
02609                 if(instanceData->objectives == NULL){
02610                         throw ErrorClass("there is no objectives object");              
02611                 }               
02612                 if(instanceData->objectives->numberOfObjectives != number){
02613                         throw ErrorClass("input number of objective not equal to number in class");             
02614                 }
02615                 if(number == 0) return true;
02616                 int i = 0;
02617                 for(i = 0; i < number; i++)instanceData->objectives->obj[i] = new Objective();
02618                 int j = 0;
02619                 if(names != NULL){
02620                         for(i = 0; i < number; i++) instanceData->objectives->obj[i]->name = names[i];                          
02621                 }       
02622                 if(maxOrMins != NULL){
02623                         for(i = 0; i < number; i++){
02624                                 if(maxOrMins[i] == "" || (maxOrMins[i].compare("max") != 0 && maxOrMins[i].compare("min") !=0)) return false;
02625                                 instanceData->objectives->obj[i]->maxOrMin = maxOrMins[i];                      
02626                         }
02627                 }
02628                 if(constants != NULL){
02629                         for(i = 0; i < number; i++) instanceData->objectives->obj[i]->constant = constants[i];                          
02630                 }
02631                 if(weights != NULL){
02632                         for(i = 0; i < number; i++) instanceData->objectives->obj[i]->weight = weights[i];                      
02633                 }
02634                 if(objectiveCoefficients != NULL){
02635                         for(i = 0; i < number; i++){
02636                                 int n = (&objectiveCoefficients[i] == NULL || objectiveCoefficients[i]->indexes == NULL)?0:objectiveCoefficients[i]->number;            
02637                                 instanceData->objectives->obj[i]->numberOfObjCoef = n;
02638                                 if(n == 0){
02639                                         instanceData->objectives->obj[i]->coef = NULL;
02640                                 }
02641                                 else{
02642                                         instanceData->objectives->obj[i]->coef = new ObjCoef*[n];
02643                                         for(j = 0; j < n; j++){
02644                                                 instanceData->objectives->obj[i]->coef[j] = new ObjCoef();
02645                                                 instanceData->objectives->obj[i]->coef[j]->idx  = objectiveCoefficients[i]->indexes[j];
02646                                                 instanceData->objectives->obj[i]->coef[j]->value = objectiveCoefficients[i]->values[j];                         
02647                                         }                       
02648                                 }                                                       
02649                         }
02650                 }
02651                 return true;
02652         }
02653         catch(const ErrorClass& eclass){
02654                 throw ErrorClass(  eclass.errormsg); 
02655         }
02656 }//setObjectives
02657 
02658 
02659 bool OSInstance::setConstraintNumber(int number){
02660         if(number < 0) return false;
02661         if(instanceData->constraints == NULL) instanceData->constraints = new Constraints();
02662         if(number == 0){
02663                 instanceData->constraints->numberOfConstraints = 0;
02664                 instanceData->constraints->con = 0;
02665                 return true;
02666         }
02667         instanceData->constraints->numberOfConstraints = number;
02668         if(instanceData->constraints->con == 0 ){
02669                 instanceData->constraints->con = new Constraint*[number];
02670         }
02671         return true;
02672 }//setConstraintNumber
02673 
02674 bool OSInstance::addConstraint(int index, string name, double lowerBound, double upperBound, double constant) {
02675         instanceData->constraints->con[ index] = new Constraint();
02676         if(index < 0 || instanceData->constraints->numberOfConstraints <= 0 || index >= instanceData->constraints->numberOfConstraints) return false;
02677         instanceData->constraints->con[ index]->name = name;
02678         if(lowerBound != -OSDBL_MAX && lowerBound != -OSDBL_MAX) instanceData->constraints->con[ index]->lb = lowerBound;
02679         if(upperBound != OSDBL_MAX && upperBound != OSDBL_MAX)instanceData->constraints->con[ index]->ub = upperBound;
02680         instanceData->constraints->con[ index]->constant = constant;
02681         return true;
02682 }//addConstraint
02683 
02684 
02685 bool OSInstance::setConstraints(int number, string* names, double* lowerBounds, double* upperBounds, double* constants){
02686         if(number < 0) return false;
02687         if(number == 0){
02688                 // this is done in setConstraintNumber
02689                 //instanceData->constraints = new Constraints();
02690                 //instanceData->constraints->numberOfConstraints = 0;
02691                 //instanceData->constraints->con = NULL;
02692                 return true;
02693         }
02694         try{
02695                 
02696                 if(instanceData->constraints  == NULL){
02697                         throw ErrorClass("there is no constraints object");             
02698                 }               
02699                 if(instanceData->constraints->numberOfConstraints != number){
02700                         throw ErrorClass("input number of constraints not equal to number in class");           
02701                 }
02702                 int i = 0;
02703                 for(i = 0; i < number; i++){
02704                         instanceData->constraints->con[i] = new Constraint();
02705                 }
02706                 if(names != NULL){
02707                         for(i = 0; i < number; i++) instanceData->constraints->con[i]->name = names[i];                         
02708                 }
02709                 if(lowerBounds != NULL){
02710                         for(i = 0; i < number; i++){
02711                                 if(lowerBounds[i] != -OSDBL_MAX && lowerBounds[i] != -OSDBL_MAX)instanceData->constraints->con[i]->lb = lowerBounds[i]; 
02712                         }
02713                 }
02714                 if(upperBounds != NULL){
02715                         for(i = 0; i < number; i++){
02716                                 if(upperBounds[i] != OSDBL_MAX && upperBounds[i] != OSDBL_MAX)instanceData->constraints->con[i]->ub = upperBounds[i]; 
02717                         }
02718                 }   
02719                 if(constants != NULL){
02720                         for(i = 0; i < number; i++) instanceData->constraints->con[i]->constant = constants[i];                         
02721                 }
02722                 return true;
02723         }
02724         catch(const ErrorClass& eclass){
02725                 throw ErrorClass(  eclass.errormsg); 
02726         }
02727 }//setConstraints
02728 
02729 bool OSInstance::setLinearConstraintCoefficients(int numberOfValues, bool isColumnMajor, 
02730                 double* values, int valuesBegin, int valuesEnd, 
02731                 int* indexes, int indexesBegin, int indexesEnd,                         
02732                 int* starts, int startsBegin, int startsEnd)
02733 {
02734         if(numberOfValues < 0) return false;
02735         if(instanceData->linearConstraintCoefficients == NULL) instanceData->linearConstraintCoefficients = new LinearConstraintCoefficients() ;
02736         if(numberOfValues == 0) return true;
02737         if((values == 0 ) ||
02738            (valuesBegin < 0 || (valuesEnd - valuesBegin + 1) != numberOfValues) ||
02739            (indexes == 0) ||
02740            (indexesBegin < 0 || (indexesEnd - indexesBegin + 1) != numberOfValues) ||
02741            (starts == 0 ) ||
02742            (startsBegin < 0  || startsBegin >= startsEnd)) return false;
02743         instanceData->linearConstraintCoefficients->numberOfValues = numberOfValues;
02744         int i = 0;
02745         //starts
02746         if(instanceData->linearConstraintCoefficients->start == NULL) instanceData->linearConstraintCoefficients->start = new IntVector();
02747         if(startsBegin == 0 ){
02748                 instanceData->linearConstraintCoefficients->start->el = starts;
02749         }
02750         else{
02751                 instanceData->linearConstraintCoefficients->start->el = new int[startsEnd - startsBegin + 1];
02752                 int k = 0;
02753                 for(i = startsBegin; i <= startsEnd; i++){
02754                         instanceData->linearConstraintCoefficients->start->el[k] = starts[i];
02755                         k++;
02756                 }
02757         }                
02758         //values
02759         if(instanceData->linearConstraintCoefficients->value == NULL) instanceData->linearConstraintCoefficients->value = new DoubleVector();
02760         if(valuesBegin == 0 ){
02761                 instanceData->linearConstraintCoefficients->value->el = values;
02762         }
02763         else{
02764                 instanceData->linearConstraintCoefficients->value->el = new double[numberOfValues];
02765                 int k = 0;
02766                 for(i = valuesBegin; i <= valuesEnd; i++){
02767                         instanceData->linearConstraintCoefficients->value->el[k] = values[i];
02768                         k++;
02769                 }
02770         }
02771         //indexes
02772         if(instanceData->linearConstraintCoefficients->rowIdx == NULL) instanceData->linearConstraintCoefficients->rowIdx = new IntVector();
02773         if(instanceData->linearConstraintCoefficients->colIdx == NULL) instanceData->linearConstraintCoefficients->colIdx = new IntVector();
02774         if(isColumnMajor){
02775                 if(indexesBegin == 0 ){
02776                         instanceData->linearConstraintCoefficients->rowIdx->el = indexes;
02777                 }
02778                 else{
02779                         instanceData->linearConstraintCoefficients->rowIdx->el = new int[numberOfValues];
02780                         int k = 0;
02781                         for(i = indexesBegin; i <= indexesEnd; i++){
02782                                 instanceData->linearConstraintCoefficients->rowIdx->el[k] = indexes[i];
02783                                 k++;
02784                         }
02785                 }
02786         } 
02787         else{
02788                 if(indexesBegin == 0 ){
02789                         instanceData->linearConstraintCoefficients->colIdx->el = indexes;
02790                 }
02791                 else{
02792                         instanceData->linearConstraintCoefficients->colIdx->el = new int[numberOfValues];
02793                         int k = 0;
02794                         for(i = indexesBegin; i <= indexesEnd; i++){
02795                                 instanceData->linearConstraintCoefficients->colIdx->el[k] = indexes[i];
02796                                 k++;
02797                         }
02798                 }
02799         }
02800         return true;
02801 }//setLinearConstraintCoefficients
02802 
02803 bool OSInstance::setQuadraticTerms(int number, 
02804                 int* rowIndexes, int* varOneIndexes, int* varTwoIndexes, double* coefficients,
02805                 int begin, int end){
02806         if(number < 0) return false;
02807         if(number != (end - begin) + 1) return false;
02808         if(number == 0){
02809                 instanceData->quadraticCoefficients = 0;
02810                 return true;
02811         }
02812         if( ((end - begin + 1) != number) ||
02813            (rowIndexes == 0) ||                    
02814            (varOneIndexes == 0) ||
02815            (varTwoIndexes == 0) ||
02816            (coefficients == 0) ) return false;
02817         instanceData->quadraticCoefficients = new QuadraticCoefficients();
02818         instanceData->quadraticCoefficients->numberOfQuadraticTerms = number;
02819         int i = 0;
02820         instanceData->quadraticCoefficients->qTerm = new QuadraticTerm*[number];
02821         for(i = 0; i < number; i++) instanceData->quadraticCoefficients->qTerm[i] = new QuadraticTerm();
02822         int k = 0;
02823         for(i = begin; i <= end; i++){
02824                         instanceData->quadraticCoefficients->qTerm[k]->idx = rowIndexes[i];
02825                         instanceData->quadraticCoefficients->qTerm[k]->idxOne = varOneIndexes[i];
02826                         instanceData->quadraticCoefficients->qTerm[k]->idxTwo = varTwoIndexes[i];
02827                         instanceData->quadraticCoefficients->qTerm[k]->coef = coefficients[i];
02828                         k++;
02829         }
02830         return true;
02831 }//setQuadraticTerms
02832 
02833 bool OSInstance::setQuadraticTermsInNonlinearExpressions(int numQPTerms, int* rowIndexes, int* varOneIndexes, int* varTwoIndexes, double* coefficients){
02834                 instanceData->nonlinearExpressions->numberOfNonlinearExpressions = numQPTerms;
02835                 instanceData->nonlinearExpressions->nl = new Nl*[ numQPTerms ];
02836                 // define the vectors
02837                 OSnLNode *nlNodePoint;
02838                 OSnLNodeVariable *nlNodeVariablePoint;
02839                 std::vector<OSnLNode*> nlNodeVec;
02840                 //
02841                 //
02842                 int i;
02843                 for(i = 0; i < numQPTerms; i++){
02844                         instanceData->nonlinearExpressions->nl[ i] = new Nl();
02845                         instanceData->nonlinearExpressions->nl[ i]->idx = rowIndexes[ i];
02846                         instanceData->nonlinearExpressions->nl[ i]->osExpressionTree = new OSExpressionTree();
02847                 // create a variable nl node for x0
02848                 nlNodeVariablePoint = new OSnLNodeVariable();
02849                 nlNodeVariablePoint->idx = varOneIndexes[ i];
02850                 // give this variable the coefficient
02851                 nlNodeVariablePoint->coef = coefficients[ i];
02852                 nlNodeVec.push_back( nlNodeVariablePoint);
02853                 // create the nl node for x1
02854                 nlNodeVariablePoint = new OSnLNodeVariable();
02855                 nlNodeVariablePoint->idx = varTwoIndexes[ i];
02856                 nlNodeVec.push_back( nlNodeVariablePoint);
02857                 // create the nl node for *
02858                 nlNodePoint = new OSnLNodeTimes();
02859                 nlNodeVec.push_back( nlNodePoint);
02860                 // the vectors are in postfix format
02861                 // now the expression tree
02862                 instanceData->nonlinearExpressions->nl[ i]->osExpressionTree->m_treeRoot =
02863                         nlNodeVec[ 0]->createExpressionTreeFromPostfix( nlNodeVec);
02864                 nlNodeVec.clear();
02865                 }
02866         return true;
02867 }//setQuadraticTermsInNonlinearExpressions
02868 
02869 bool OSInstance::initializeNonLinearStructures( ){
02870         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02871         if( m_bNonLinearStructuresInitialized == true) return true;
02872         if( m_bProcessVariables == false) processVariables();
02873         if( m_bProcessObjectives == false) processObjectives();
02874         if( m_bProcessConstraints == false) processConstraints();
02875         m_iVariableNumber = instanceData->variables->numberOfVariables;
02876         m_iConstraintNumber = instanceData->constraints->numberOfConstraints;
02877         m_iObjectiveNumber = instanceData->objectives->numberOfObjectives;
02878         // get all of the expression trees
02879         if( m_bProcessExpressionTrees == false) getAllNonlinearExpressionTrees();
02880         // before proceeding get a copy of the map of the Expression Trees
02881         if( m_bDuplicateExpressionTreesMap == false) duplicateExpressionTreesMap();
02882         // now create all of the variable maps for each expression tree
02883         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){  
02884                 (posMapExpTree->second)->getVariableIndiciesMap() ;
02885         }
02886         // add the quadratic terms if necessary
02887         if(getNumberOfQuadraticTerms() > 0) addQTermsToExressionTree();
02888         // now get the map of all nonlinear variables
02889         getAllNonlinearVariablesIndexMap( );
02890         getDenseObjectiveCoefficients();
02891         m_mdConstraintFunctionValues = new double[ this->instanceData->constraints->numberOfConstraints];
02892         m_mdObjectiveFunctionValues = new double[ this->instanceData->objectives->numberOfObjectives];
02893         //m_mdObjGradient = new double[ this->instanceData->variables->numberOfVariables];
02894         m_bNonLinearStructuresInitialized = true;
02895         m_bProcessVariables = true;
02896         m_bProcessObjectives = true;
02897         m_bProcessConstraints = true;
02898         m_bProcessExpressionTrees = true;
02899         m_bDuplicateExpressionTreesMap = true;
02900         return true;
02901 }
02902 
02903 SparseJacobianMatrix *OSInstance::getJacobianSparsityPattern( ){
02904         // if already called return the sparse Jacobian
02905         // it is important that this method NOT get called twice -- if
02906         // there are linear terms in <linearConstraintCoefficients> that
02907         // also appear in <nonlinearExpressions> then they will keep getting added
02908         // to the modified expession tree with each call to this method
02909         if( m_bSparseJacobianCalculated == true) return m_sparseJacMatrix;
02910         //std::cout << "INSIDE GET JACOBIAN SPARSITY PATTERN" << std::endl;
02911         // determine if we are in column or row major
02912         getLinearConstraintCoefficientMajor();
02913         // make sure the data structures have been inialized
02914         if( m_bNonLinearStructuresInitialized == false) initializeNonLinearStructures( );
02915         try{
02916                 if( m_bColumnMajor == true){
02917                          if( getSparseJacobianFromColumnMajor( ) == false) throw ErrorClass("An error occurred in getSpareJacobianFromColumnMajor");
02918                 }
02919                 else {
02920                         if( getSparseJacobianFromRowMajor( ) == false) throw ErrorClass("An error occurred in getSpareJacobianFromRowMajor");
02921                 }
02922         }
02923         catch(const ErrorClass& eclass){
02924                 throw ErrorClass(  eclass.errormsg); 
02925         }
02926         // now fill in the arrays of the sparseJacMatrix
02927         m_sparseJacMatrix = new SparseJacobianMatrix();
02928         // we point to memory already created so don't 
02929         // destroy during garbage collection
02930         m_sparseJacMatrix->bDeleteArrays = false;
02931         m_sparseJacMatrix->valueSize =  m_iJacValueSize;
02932         m_sparseJacMatrix->starts = m_miJacStart;
02933         m_sparseJacMatrix->conVals = m_miJacNumConTerms;
02934         m_sparseJacMatrix->indexes = m_miJacIndex;
02935         m_sparseJacMatrix->values = m_mdJacValue;
02936         m_bSparseJacobianCalculated = true;
02937         return m_sparseJacMatrix;
02938 }//getJacobianSparsityPatter
02939 
02940 bool OSInstance::addQTermsToExressionTree(){
02941         int i, k, idx;
02942         // get the number of qTerms
02943         int numQTerms = instanceData->quadraticCoefficients->numberOfQuadraticTerms;    
02944         if(numQTerms <= 0 || m_bQTermsAdded == true) return true;
02945         OSnLNodeVariable* nlNodeVariableOne;
02946         OSnLNodeVariable* nlNodeVariableTwo;
02947         OSnLNodeTimes* nlNodeTimes;
02948         OSnLNodePlus* nlNodePlus;
02949         OSExpressionTree* expTree;
02950         getQuadraticTerms();    
02951         //std::cout << "PROCESSING QUADRATIC TERMS" << std::endl;
02952         for(i = 0; i < numQTerms; i++){
02953                 idx = m_quadraticTerms->rowIndexes[ i];
02954                 //std::cout << "PROCESSING QTERM  = "  << i <<std::endl;
02955                 // see if row idx is in the expression tree
02956                 if( m_mapExpressionTreesMod.find( idx) != m_mapExpressionTreesMod.end() ) {
02957                         // row idx is in the expression tree
02958                         // add the qTerm in row idx  to the expression tree     
02959                         // define two new OSnLVariable nodes, an OSnLnodeTimes, and OSnLnodePlus 
02960                         nlNodeVariableOne = new OSnLNodeVariable();
02961                         nlNodeVariableOne->idx = m_quadraticTerms->varOneIndexes[ i];
02962                         // see if the variable indexed by nlNodeVariableOne->idx is in the expression tree for row idx
02963                         // if not, add to mapVarIdx
02964                         expTree = m_mapExpressionTreesMod[ idx];
02965                         if(  expTree->m_bIndexMapGenerated == false) expTree->getVariableIndiciesMap();
02966                         if( (*expTree->mapVarIdx).find( nlNodeVariableOne->idx) == (*expTree->mapVarIdx).end()  ){
02967                                 // add to map
02968                                 k = (*expTree->mapVarIdx).size();
02969                                 (*expTree->mapVarIdx)[ nlNodeVariableOne->idx] =  k + 1;
02970                                 #ifdef DEBUG
02971                                 std::cout << "ADDED THE FOLLOWING VAIRABLE TO THE MAP" <<  nlNodeVariableOne->idx << std::endl;
02972                                 #endif
02973                         }
02974                         nlNodeVariableOne->coef = m_quadraticTerms->coefficients[ i];
02975                         nlNodeVariableTwo = new OSnLNodeVariable();
02976                         nlNodeVariableTwo->idx = m_quadraticTerms->varTwoIndexes[ i];
02977                         // see if the variable indexed by nlNodeVariableTwo->idx is in the expression tree for row idx
02978                         // if not, add to mapVarIdx
02979                         if( (*expTree->mapVarIdx).find( nlNodeVariableTwo->idx) == (*expTree->mapVarIdx).end()  ){
02980                                 // add to map
02981                                 k = (*expTree->mapVarIdx).size();
02982                                 (*expTree->mapVarIdx)[ nlNodeVariableTwo->idx] =  k + 1;
02983                                 #ifdef DEBUG
02984                                 std::cout << "ADDED THE FOLLOWING VAIRABLE TO THE MAP" <<  nlNodeVariableTwo->idx << std::endl;
02985                                 #endif
02986                         }
02987                         nlNodeVariableTwo->coef = 1.;
02988                         // now multiply the two new variable nodes together
02989                         nlNodeTimes = new OSnLNodeTimes();
02990                         nlNodeTimes->m_mChildren[ 0] = nlNodeVariableOne;
02991                         nlNodeTimes->m_mChildren[ 1] = nlNodeVariableTwo;               
02992                         // now add the result to the expression tree
02993                         nlNodePlus = new OSnLNodePlus();
02994                         nlNodePlus->m_mChildren[ 0] = expTree->m_treeRoot;
02995                         nlNodePlus->m_mChildren[ 1] = nlNodeTimes;
02996                         //expTree = new OSExpressionTree();
02997                         expTree->m_treeRoot = nlNodePlus ;
02998                         // get rid of old variable map
02999                         if(expTree->m_bIndexMapGenerated == true){
03000                                 delete expTree->mapVarIdx;
03001                                 expTree->mapVarIdx = NULL;
03002                                 expTree->m_bIndexMapGenerated = false;
03003                         }       
03004                         //expTree->m_bIndexMapGenerated = false;
03005                         //m_mapExpressionTreesMod[ idx ]  = expTree;    
03006                         //expTree->mapVarIdx = m_mapExpressionTreesMod[ idx]->mapVarIdx;
03007                 }
03008                 else{
03009                         // create the quadratic expression to add to the expression tree
03010                         nlNodeVariableOne = new OSnLNodeVariable();
03011                         nlNodeVariableOne->idx = m_quadraticTerms->varOneIndexes[ i];
03012                         nlNodeVariableOne->coef = m_quadraticTerms->coefficients[ i];
03013                         nlNodeVariableTwo = new OSnLNodeVariable();
03014                         nlNodeVariableTwo->idx = m_quadraticTerms->varTwoIndexes[ i];
03015                         nlNodeVariableTwo->coef = 1.;
03016                         // now multiply the two new variable nodes together
03017                         nlNodeTimes = new OSnLNodeTimes();
03018                         nlNodeTimes->m_mChildren[ 0] = nlNodeVariableOne;
03019                         nlNodeTimes->m_mChildren[ 1] = nlNodeVariableTwo;
03020                         // create a new expression tree corresponding to row idx.
03021                         expTree = new OSExpressionTree();                                               
03022                         expTree->m_treeRoot = nlNodeTimes ;
03023                         expTree->mapVarIdx = expTree->getVariableIndiciesMap();         
03024                         m_mapExpressionTreesMod[ idx ]  = expTree;
03025                         if(idx < 0){
03026                                 m_iObjectiveNumberNonlinear++;
03027                                 m_bProcessExpressionTrees = true;
03028                         }
03029                         else{
03030                                 m_iConstraintNumberNonlinear++;
03031                                 m_bProcessExpressionTrees = true;
03032                         }
03033                         #ifdef DEBUG
03034                         std::cout << "NUMBER OF EXPRESSION TREES = "  << m_mapExpressionTreesMod.size() <<std::endl;
03035                         std::cout << "NUMBER OF NONLINEAR OBJECTIVES = "  << getNumberOfNonlinearObjectives() <<std::endl;
03036                         #endif
03037                 } 
03038                 // if there were no nonlinear terms make this the expression tree
03039                 if(m_iNonlinearExpressionNumber <= 0) m_mapExpressionTrees = m_mapExpressionTreesMod;
03040                 m_bQTermsAdded =true;
03041         }
03042         return true;
03043 }
03044 
03045 double OSInstance::calculateFunctionValue(int idx, double *x, bool new_x){
03046         try{
03047 
03048                 int i, j;
03049                 double dvalue = 0;
03050                 if( m_binitForAlgDiff == false) initForAlgDiff();
03051                 if( m_bSparseJacobianCalculated == false) getJacobianSparsityPattern();
03052                 if(idx >= 0){ // we have a constraint
03053                         // make sure the index idx is valid
03054                         if( getConstraintNumber() <= idx  ) throw 
03055                         ErrorClass("constraint index not valid in OSInstance::calculateFunctionValue");
03056                         if( new_x == false) return *(m_mdConstraintFunctionValues + idx);
03057                         // get the nonlinear part
03058                         if( m_mapExpressionTreesMod.find( idx) != m_mapExpressionTreesMod.end() ){
03059                                 dvalue = m_mapExpressionTreesMod[ idx]->calculateFunction( x,  new_x);
03060                                 //dvalue = vdFunVals[ idx + 1];
03061                         }
03062                         // now the linear part
03063                         // be careful, loop over only the constant terms in sparseJacMatrix
03064                         i = m_sparseJacMatrix->starts[ idx];
03065                         j = m_sparseJacMatrix->starts[ idx + 1 ];
03066                         //while ( i <  j &&  (i - m_sparseJacMatrix->starts[ idx])  < m_sparseJacMatrix->conVals[ idx] ){
03067                         while ( (i - m_sparseJacMatrix->starts[ idx])  < m_sparseJacMatrix->conVals[ idx] ){
03068                                 //std::cout << "m_sparseJacMatrix->values[ i] " << m_sparseJacMatrix->values[ i] << std::endl;
03069                                 //std::cout << "m_sparseJacMatrix->indexes[ i] " << m_sparseJacMatrix->indexes[ i] << std::endl;
03070                                 dvalue += m_sparseJacMatrix->values[ i]*x[ m_sparseJacMatrix->indexes[ i] ];
03071                                 i++;
03072                         }       
03073                         // add in the constraint function constant
03074                         dvalue += m_mdConstraintConstants[ idx ];
03075                         return dvalue;
03076                 }
03077                 else{ // we have an objective function
03078                         // make sure the index idx is valid
03079                         if( getObjectiveNumber() <= ( abs( idx) - 1) ) throw 
03080                         ErrorClass("objective function index not valid in OSInstance::calculateFunctionValue");
03081                         if( new_x == false) return *(m_mdObjectiveFunctionValues + ( abs( idx) - 1));
03082                         // get the nonlinear part
03083                         if( m_mapExpressionTreesMod.find( idx) != m_mapExpressionTreesMod.end() ){
03084                                 dvalue = m_mapExpressionTreesMod[ idx]->calculateFunction( x,  new_x);
03085                         }
03086                         // get linear part
03087                         SparseVector **objCoef = getObjectiveCoefficients();
03088                         SparseVector *obj = objCoef[ abs( idx) - 1];
03089                         for(i = 0; i < obj->number; i++){
03090                                 dvalue += x[ obj->indexes[i]]*(obj->values[ i]);
03091                         }
03092                         // add in the objective function constant
03093                         dvalue += m_mdObjectiveConstants[ abs( idx) - 1 ];
03094                         // get the coefficients for objective function idx
03095                         *(m_mdObjectiveFunctionValues + ( abs( idx) - 1)) = dvalue;
03096                         return *(m_mdObjectiveFunctionValues + ( abs( idx) - 1));
03097                 }
03098         }
03099         catch(const ErrorClass& eclass){
03100                 throw ErrorClass( eclass.errormsg);
03101         } 
03102 }//calculateFunctionValue
03103 
03104 
03105 double *OSInstance::calculateAllConstraintFunctionValues( double* x, double *objLambda, double *conLambda,
03106         bool new_x, int highestOrder){  
03107         try{
03108                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) 
03109                         getIterateResults(x, objLambda, conLambda, new_x,  highestOrder);
03110         }
03111         catch(const ErrorClass& eclass){
03112                 throw ErrorClass( eclass.errormsg);
03113         } 
03114         return m_mdConstraintFunctionValues;
03115 }//calculateAllConstraintFunctionValues
03116 
03117 
03118 double *OSInstance::calculateAllConstraintFunctionValues(double* x, bool new_x){
03119         try{
03120                 m_iHighestOrderEvaluated = -1;
03121                 if( new_x == false) return m_mdConstraintFunctionValues;
03122                 int idx, numConstraints;
03123                 numConstraints = getConstraintNumber();
03124                 // loop over all constraints
03125                 for(idx = 0; idx < numConstraints; idx++){
03126                         m_mdConstraintFunctionValues[ idx]  = calculateFunctionValue(idx, x, new_x);    
03127                 }
03128                 
03129         }
03130         catch(const ErrorClass& eclass){
03131                 throw ErrorClass( eclass.errormsg);
03132         } 
03133         return m_mdConstraintFunctionValues;    
03134 }//end calculateAllConstraintFunctionValues
03135 
03136 
03137 double *OSInstance::calculateAllObjectiveFunctionValues( double* x, double *objLambda, double *conLambda,
03138         bool new_x, int highestOrder){  
03139         try{
03140                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) 
03141                         getIterateResults(x, objLambda, conLambda, new_x,  highestOrder);
03142         }
03143         catch(const ErrorClass& eclass){
03144                 throw ErrorClass( eclass.errormsg);
03145         } 
03146         return m_mdObjectiveFunctionValues;
03147 }//calculateAllConstraintFunctionValues
03148 
03149 
03150 double *OSInstance::calculateAllObjectiveFunctionValues( double* x, bool new_x){        
03151         try{
03152                 m_iHighestOrderEvaluated = -1;
03153                 
03154                 if( new_x == false) return m_mdObjectiveFunctionValues;
03155                 int idx, numObjectives;
03156                 numObjectives = getObjectiveNumber();
03157                 // loop over all constraints
03158                 for(idx = 0; idx < numObjectives; idx++){
03159                         m_mdObjectiveFunctionValues[ idx]  = calculateFunctionValue(-idx -1, x, new_x); 
03160                 }
03161         }
03162         catch(const ErrorClass& eclass){
03163                 throw ErrorClass( eclass.errormsg);
03164         } 
03165         return m_mdObjectiveFunctionValues;
03166 }//calculateAllObjectiveFunctionValues
03167 
03168 
03169 SparseJacobianMatrix *OSInstance::calculateAllConstraintFunctionGradients(double* x, double *objLambda, double *conLambda,
03170                 bool new_x, int highestOrder){
03171         try{
03172                 if(highestOrder < 1 ) throw ErrorClass("When calling calculateAllConstraintFunctionGradients highestOrder should be 1 or 2");
03173                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) 
03174                         getIterateResults(x, objLambda, conLambda,  new_x,  highestOrder);
03175         }//end try
03176         catch(const ErrorClass& eclass){
03177                 throw ErrorClass( eclass.errormsg);
03178         } 
03179         return m_sparseJacMatrix;
03180 }//calculateAllConstraintFunctionGradients      
03181 
03182 
03183 
03184 SparseVector *OSInstance::calculateConstraintFunctionGradient(double* x, double *objLambda, double *conLambda,
03185                 int idx, bool new_x, int highestOrder){
03186         try{
03187                 if(highestOrder < 1 ) throw ErrorClass("When calling calculateConstraintFunctionGradient highestOrder should be 1 or 2");
03188                 if(idx < 0 || idx >= instanceData->constraints->numberOfConstraints ) 
03189                         throw ErrorClass("invalid index passed to calculateConstraintFunctionGrad");
03190                 SparseVector *sp;
03191                 sp = new SparseVector();
03192                 sp->bDeleteArrays = true;
03193                 int i;
03194                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) 
03195                         getIterateResults(x, objLambda, conLambda,  new_x,  highestOrder);
03196                 sp->number = m_miJacStart[ idx + 1] - m_miJacStart[ idx];
03197                 sp->values = new double[ sp->number];
03198                 sp->indexes = new int[ sp->number];
03199                 for(i = 0; i < sp->number; i++){
03200                         sp->values[ i] = m_mdJacValue[ m_miJacStart[ idx] +  i];
03201                         sp->indexes[ i] = m_miJacIndex[ m_miJacStart[ idx] +  i];
03202                 }
03203                 return sp;
03204         }
03205         catch(const ErrorClass& eclass){
03206                 throw ErrorClass( eclass.errormsg);
03207         } 
03208 }//calculateConstraintFunctionGradient
03209 
03210 
03211 SparseVector *OSInstance::calculateConstraintFunctionGradient(double* x, int idx, bool new_x){
03212         try{
03213                 if(idx < 0 || idx >= instanceData->constraints->numberOfConstraints ) 
03214                         throw ErrorClass("invalid index passed to calculateConstraintFunctionGrad");
03215                 SparseVector *sp;
03216                 sp = new SparseVector();
03217                 sp->bDeleteArrays = true;
03218                 int i;
03219                 if( new_x == true || (1 > m_iHighestOrderEvaluated)  ) 
03220                         getIterateResults(x, NULL, NULL,  new_x,  1);
03221                 sp->number = m_miJacStart[ idx + 1] - m_miJacStart[ idx];
03222                 sp->values = new double[ sp->number];
03223                 sp->indexes = new int[ sp->number];
03224                 for(i = 0; i < sp->number; i++){
03225                         sp->values[ i] = m_mdJacValue[ m_miJacStart[ idx] +  i];
03226                         sp->indexes[ i] = m_miJacIndex[ m_miJacStart[ idx] + i];
03227                 }
03228                 return sp;
03229         }
03230         catch(const ErrorClass& eclass){
03231                 throw ErrorClass( eclass.errormsg);
03232         } 
03233 }//calculateConstraintFunctionGradient
03234 
03235 
03236 double **OSInstance::calculateAllObjectiveFunctionGradients(double* x, double *objLambda, double *conLambda,
03237                 bool new_x, int highestOrder){
03238         try{
03239                 if(highestOrder < 1 ) throw ErrorClass("When calling calculateAllObjectiveFunctionGradients highestOrder should be 1 or 2");
03240                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) {
03241                         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03242                         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03243                                 if(posMapExpTree->first < 0){ // this nonlinear expression indexes an objective function
03244                                         m_mmdObjGradient[ abs( posMapExpTree->first) - 1 ] = calculateObjectiveFunctionGradient(x, objLambda, conLambda, 
03245                                                         posMapExpTree->first, new_x, highestOrder);
03246                                 }       
03247                         }
03248                 }
03249         }
03250         catch(const ErrorClass& eclass){
03251                 throw ErrorClass( eclass.errormsg);
03252         } 
03253         return m_mmdObjGradient;
03254 }// calculateAllObjectiveFunctionGradients
03255 
03256 double *OSInstance::calculateObjectiveFunctionGradient(double* x, double *objLambda, double *conLambda,
03257                 int objIdx, bool new_x, int highestOrder){
03258         /* if we are just doing an objective function gradient we should do a zero order forward sweep
03259          * and a reverse first order reverse sweep
03260          */
03261         try{
03262                 if(highestOrder < 1 ) throw ErrorClass("When calling calculateObjectiveFunctionGradient highestOrder should be 1 or 2");        
03263                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) {
03264                         int domainIdx = 0;      
03265                         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03266                         std::map<int, int>::iterator posVarIndexMap;            
03267                         int iHighestOrderEvaluatedStore;
03268                         unsigned int i;
03269                         iHighestOrderEvaluatedStore = m_iHighestOrderEvaluated;
03270                         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03271                                 //kipp: modify for more than one obj
03272                                 if(posMapExpTree->first == objIdx){
03273                                         if( new_x == true ){
03274                                                 if( m_vdX.size() > 0) m_vdX.clear();
03275                                                 for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
03276                                                         m_vdX.push_back( x[ posVarIndexMap->first]) ;
03277                                                 }
03278                                                 if( (m_bOSADFunIsCreated == false || m_bCppADMustReTape == true )  && (m_mapExpressionTreesMod.size() > 0) ) {
03279                                                         createOSADFun( m_vdX);
03280                                                 }                                               
03281                                         }       
03282                                         
03283                                         if(( new_x == true ) || (m_iHighestOrderEvaluated < 0) )this->forwardAD(0, m_vdX);
03284                                         
03285                                         if(( new_x == true ) || (m_iHighestOrderEvaluated < 1) ) {
03286                                                 m_vdRangeUnitVec[ domainIdx] = 1.;
03287                                                 m_vdYjacval = this->reverseAD(1, m_vdRangeUnitVec);
03288                                                 for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
03289                                                         m_mmdObjGradient[ abs( objIdx) - 1 ][ m_miNonLinearVarsReverseMap[ i]] = m_vdYjacval[ i] + 
03290                                                                 m_mmdDenseObjectiveCoefficients[  abs( objIdx) - 1][ m_miNonLinearVarsReverseMap[ i]];
03291                                                 }
03292                                         }
03293                                         m_iHighestOrderEvaluated = iHighestOrderEvaluatedStore;
03294                                         m_vdRangeUnitVec[ domainIdx] = 0.;
03295                                         // exit the loop
03296                                         break;
03297                                 }
03298                                 domainIdx++;
03299                         }
03300                 }                       
03301         }
03302         catch(const ErrorClass& eclass){
03303                 throw ErrorClass( eclass.errormsg);
03304         } 
03305         return m_mmdObjGradient[abs( objIdx) - 1];
03306 }//calculateObjectiveFunctionGradient
03307 
03308 
03309 double *OSInstance::calculateObjectiveFunctionGradient(double* x, int objIdx, bool new_x){
03310         try{
03311                 int domainIdx = 0;      
03312                 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03313                 std::map<int, int>::iterator posVarIndexMap;    
03314                 unsigned int i;
03315                 int  iHighestOrderEvaluatedStore;       
03316                 iHighestOrderEvaluatedStore = m_iHighestOrderEvaluated;
03317                 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03318                         if(posMapExpTree->first == objIdx){
03319                                 if( new_x == true ){
03320                                         if( m_vdX.size() > 0) m_vdX.clear();
03321                                         for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
03322                                                 m_vdX.push_back( x[ posVarIndexMap->first]) ;
03323                                         }
03324                                         if( (m_bOSADFunIsCreated == false || m_bCppADMustReTape == true )  && (m_mapExpressionTreesMod.size() > 0) ) {
03325                                                 createOSADFun( m_vdX);
03326                                         }                                               
03327                                 }       
03328                                 
03329                                 if(( new_x == true ) || (m_iHighestOrderEvaluated < 0) )this->forwardAD(0, m_vdX);
03330                                 
03331                                 if(( new_x == true ) || (m_iHighestOrderEvaluated < 1) ) {
03332                                         m_vdRangeUnitVec[ domainIdx] = 1.;
03333                                         m_vdYjacval = this->reverseAD(1, m_vdRangeUnitVec);
03334                                         for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
03335                                                 m_mmdObjGradient[ abs( objIdx) - 1 ][ m_miNonLinearVarsReverseMap[ i]] = m_vdYjacval[ i] + 
03336                                                         m_mmdDenseObjectiveCoefficients[  abs( objIdx) - 1][ m_miNonLinearVarsReverseMap[ i]];
03337                                         }
03338                                 }
03339                                 m_iHighestOrderEvaluated = iHighestOrderEvaluatedStore;
03340                                 m_vdRangeUnitVec[ domainIdx] = 0.;
03341                                 // exit the loop
03342                                 break;
03343                         }
03344                         domainIdx++;
03345                 }
03346         }       
03347         catch(const ErrorClass& eclass){
03348                 throw ErrorClass( eclass.errormsg);
03349         } 
03350         return m_mmdObjGradient[abs( objIdx) - 1];
03351 }//calculateObjectiveFunctionGradient   
03352 
03353 SparseHessianMatrix *OSInstance::calculateLagrangianHessian( double* x, double *objLambda, double *conLambda,
03354         bool new_x, int highestOrder){
03355         try{
03356                 if(highestOrder != 2 ) throw ErrorClass("When calling calculateLagrangianHessian highestOrder should be 2");
03357                 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated)  ) {
03358                         //std::cout  << "CALL getIterateResults() FROM calculateLagrangianHessain" << std::endl;
03359                         getIterateResults(x, objLambda, conLambda,  new_x,  highestOrder);
03360                 }
03361         }
03362         catch(const ErrorClass& eclass){
03363                 throw ErrorClass( eclass.errormsg);
03364         } 
03365         return m_LagrangianSparseHessian;
03366 }//calculateLagrangianHessian
03367 
03368 SparseHessianMatrix *OSInstance::calculateHessian(double* x, int idx, bool new_x){
03369         try{
03370                 if(idx < -instanceData->objectives->numberOfObjectives || idx >= instanceData->constraints->numberOfConstraints ) 
03371                         throw ErrorClass("invalid index passed to calculateHessian");
03372                 double *objLambda = new double[ getObjectiveNumber() ];
03373                 double *conLambda = new double[ getConstraintNumber() ];
03374                 //std::cout << "NUMBER OF OBJECTIVES = " << getObjectiveNumber() << std::endl;
03375                 //std::cout << "NUMBER OF CONSTRAINTS = " << getConstraintNumber() << std::endl;
03376                 int i;
03377                 // initialize all to zero
03378                 for(i = 0; i < getObjectiveNumber(); i++){
03379                         objLambda[ i] = 0.0;
03380                 }
03381                 for(i = 0; i < getConstraintNumber(); i++){
03382                         conLambda[ i] = 0.0;
03383                 }
03384                 // see if we have the index of an objective function or a constraint
03385                 // and more corresponding component 1.0
03386                 if(idx < 0){
03387                         objLambda[ abs(idx) - 1] = 1.0;
03388                 }
03389                 else{
03390                         conLambda[ idx] = 1.0;
03391                 }
03392                 if( new_x == true || (2 > m_iHighestOrderEvaluated)  ) {
03393                         getIterateResults(x, objLambda, conLambda, new_x,  2);
03394 //                      std::cout  << "CALL getIterateResults() FROM calculateHessian" << std::endl;
03395                 }
03396                 delete[] objLambda;
03397                 delete[] conLambda;
03398         }
03399         catch(const ErrorClass& eclass){
03400                 throw ErrorClass( eclass.errormsg);
03401         } 
03402         return m_LagrangianSparseHessian;
03403 }//calculateHessian                     
03404 
03405 
03406 
03407 bool OSInstance::getSparseJacobianFromColumnMajor( ){
03408         // we assume column major matrix
03409         if( m_bColumnMajor == false) return false;
03410         int iNumRowStarts = getConstraintNumber() + 1;  
03411         int i,j, iTemp;
03412         int iNumVariableStarts = getVariableNumber() ;
03413         int *start = NULL;
03414         int *index = NULL;
03415         double *value = NULL;
03416         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
03417                 start = this->instanceData->linearConstraintCoefficients->start->el;
03418                 index = this->instanceData->linearConstraintCoefficients->rowIdx->el;
03419                 value = this->instanceData->linearConstraintCoefficients->value->el;
03420         }
03421         m_miJacStart = new int[ iNumRowStarts];
03422         m_miJacNumConTerms = new int[ getConstraintNumber()];
03423         OSnLNodePlus *nlNodePlus;
03424         OSnLNodeVariable *nlNodeVariable;
03425         OSExpressionTree *expTree = NULL;
03426         // now initialize starts and variable index maps 
03427         for ( i = 0; i < iNumRowStarts; i++){                   
03428                 m_miJacStart [ i ] = 0;
03429                 // map the variables  in the nonlinear rows
03430                 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
03431                         // the following is equivalent to  m_treeRoot->getVariableIndexMap( i);
03432                         m_mapExpressionTreesMod[ i]->getVariableIndiciesMap();
03433                         
03434                 }
03435         }
03436         // only execute the following code if there are linear constraint coefficients
03437         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
03438                 // i is indexing columns (variables) and j is indexing row numbers 
03439                 for (i = 0; i < iNumVariableStarts; i++){       
03440                         for (j = start[i]; j < start[ i + 1 ]; j++){
03441                                 // index[ j] is a row index, we have just found an occurance of row index[j]
03442                                 // therefore we increase by 1 (or push back) the start of the row indexed by index[j] + 1, 
03443                                 // i.e. the start of the next row
03444                                 // check to see if variable i is linear/constant in the row index[ j] 
03445                                 // if so, increment m_miJacStart[ index[j] + 1]
03446                                 //
03447                                 if( (m_mapExpressionTreesMod.find( index[ j]) != m_mapExpressionTreesMod.end() ) &&
03448                                         ( (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).find( i) != (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).end()) ){
03449                                         // variable i is appears in the expression tree for row index[ j]
03450                                         // add the coefficient corresponding to variable i in row index[ j] to the expression tree      
03451                                         // define a new OSnLVariable and OSnLnodePlus 
03452                                         // don't add a zero
03453                                         if( value[j] > 0 || value[j] < 0){
03454                                                 expTree = m_mapExpressionTreesMod[ index[j]  ];
03455                                                 nlNodeVariable = new OSnLNodeVariable();
03456                                                 nlNodeVariable->coef = value[ j];
03457                                                 nlNodeVariable->idx = i;
03458                                                 nlNodePlus = new OSnLNodePlus();
03459                                                 nlNodePlus->m_mChildren[ 0] = m_mapExpressionTreesMod[ index[ j] ]->m_treeRoot;
03460                                                 nlNodePlus->m_mChildren[ 1] = nlNodeVariable;
03461                                                 //expTree = new OSExpressionTree();
03462                                                 expTree->m_treeRoot = nlNodePlus ;
03463                                                 //expTree->mapVarIdx = m_mapExpressionTreesMod[ index[ j]]->mapVarIdx;
03464                                                 //m_mapExpressionTreesMod[ index[ j] ]  = expTree;      
03465                                                 //std::cout << m_mapExpressionTreesMod[ index[ j] ]->m_treeRoot->getNonlinearExpressionInXML() << std::endl;    
03466                                                 //std::cout << m_mapExpressionTrees[ index[ j] ]->m_treeRoot->getNonlinearExpressionInXML() << std::endl;
03467                                         }
03468                                 }
03469                                 else{ 
03470                                         m_miJacStart[ index[j] + 1] ++;
03471                                 }                               
03472                         }
03473                 }
03474         }
03475         // at this point, m_miJacStart[ i] holds the number of columns with a linear/constant nonzero in row i - 1
03476         // we are not done with the start indicies, if we are here, and we
03477         // knew the correct starting point of row i -1, the correct starting point
03478         // for row i is m_miJacStart[i] + m_miJacStart [i - 1]
03479         m_miJacStart[0] = 0;
03480         for (i = 1; i < iNumRowStarts; i++ ){
03481                 m_miJacNumConTerms[ i - 1] = m_miJacStart[i];
03482                 if( m_mapExpressionTreesMod.find( i - 1) != m_mapExpressionTreesMod.end() ){
03483                         m_miJacStart[i] += (m_miJacStart[i - 1] + (*m_mapExpressionTreesMod[ i - 1]->mapVarIdx).size() );
03484                 }
03485                 else{
03486                         m_miJacStart[i] += m_miJacStart[i - 1];
03487                 }       
03488         }
03489         // dimension miIndex and mdValue here
03490         m_iJacValueSize =       m_miJacStart[ iNumRowStarts - 1];
03491         m_miJacIndex = new int[  m_iJacValueSize];
03492         m_mdJacValue = new double[ m_iJacValueSize ];
03493         // now get the values of the constant terms if there are any
03494         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
03495                 // loop over variables  
03496                 for (i = 0; i < iNumVariableStarts; i++){
03497                         // get row indices and values of the A matrix
03498                         // kipp -- should we have a check to see if start[i+1] > start[i]
03499                         for (j = start[i]; j < start[ i + 1 ]; j++){
03500                                 // store this variable index in every row where the variable appears
03501                                 // however, don't store this as constant term if it appears in mapVarIdx
03502                                 if( (m_mapExpressionTreesMod.find( index[ j]) == m_mapExpressionTreesMod.end() ) ||
03503                                         ( (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).find( i) == (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).end()) ){
03504                                         iTemp = m_miJacStart[ index[j]];
03505                                         m_miJacIndex[ iTemp] = i;
03506                                         m_mdJacValue[ iTemp] = value[j];
03507                                         m_miJacStart[ index[j]]++;      
03508                                 }               
03509                         }                       
03510                 } 
03511         }
03512         //
03513         std::map<int, int>::iterator posVarIdx;
03514         // m_miJacStart[ i] is now equal to the correct m_miJacStart[ i] + m_miJacNumConTerms[ i], so readjust
03515         for (i = 0; i < iNumRowStarts - 1; i++ ){
03516                 m_miJacStart[ i] = m_miJacStart [ i] - m_miJacNumConTerms[ i] ; 
03517                 iTemp = m_miJacStart[ i] + m_miJacNumConTerms[ i];
03518                 // if the row is in the list of expression trees read in idicies and values
03519                 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
03520                         for(posVarIdx = (*m_mapExpressionTreesMod[ i]->mapVarIdx).begin(); posVarIdx 
03521                         != (*m_mapExpressionTreesMod[ i]->mapVarIdx).end(); ++posVarIdx){
03522                                 m_miJacIndex[ iTemp] = posVarIdx->first;
03523                                 m_mdJacValue[ iTemp] = 0;
03524                                 iTemp++;
03525                         }
03526                 }
03527         }
03528         #ifdef DEBUG
03529         std::cout << "HERE ARE ROW STARTS:" << std::endl;
03530         for (i = 0; i < iNumRowStarts; i++ ){
03531                 std::cout <<  m_miJacStart[ i] << "  "; 
03532         }
03533         std::cout << std::endl << std::endl;
03534         std::cout << "HERE ARE VARIABLE INDICIES:" << std::endl;
03535         for (i = 0; i < m_miJacStart[ iNumRowStarts - 1]; i++ ){
03536                 std::cout <<  m_miJacIndex[ i] << "  "; 
03537         }
03538         std::cout << std::endl << std::endl;
03539         std::cout << "HERE ARE VALUES:" << std::endl;
03540         for (i = 0; i < m_miJacStart[ iNumRowStarts - 1]; i++ ){
03541                 std::cout <<  m_mdJacValue[ i] << "  "; 
03542         }
03543         std::cout << std::endl << std::endl;
03544 
03545         std::cout << "HERE ARE NUMBER OF CONSTANT TERMS:" << std::endl;
03546         for (i = 0; i < iNumRowStarts - 1; i++ ){
03547                 std::cout <<  m_miJacNumConTerms[ i ] << "  ";  
03548         }
03549         std::cout << std::endl << std::endl;
03550         #endif
03551         return true;
03552 }//getSparseJacobianFromColumnMajor
03553 
03554 
03555 
03556 bool OSInstance::getSparseJacobianFromRowMajor( ){
03557         // we assume row major matrix
03558         if( m_bColumnMajor == true) return false;
03559         int iNumJacRowStarts = getConstraintNumber() + 1;
03560         std::map<int, int>::iterator posVarIdx; 
03561         int i,j, k;
03562         int *start = NULL;
03563         int *index = NULL;
03564         double *value = NULL;
03565         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
03566                 start = this->instanceData->linearConstraintCoefficients->start->el;
03567                 index = this->instanceData->linearConstraintCoefficients->colIdx->el;
03568                 value = this->instanceData->linearConstraintCoefficients->value->el;
03569         }
03570         m_miJacStart = new int[ iNumJacRowStarts];
03571         m_miJacNumConTerms = new int[ getConstraintNumber()];
03572         OSnLNodePlus *nlNodePlus;
03573         OSnLNodeVariable *nlNodeVariable;
03574         //OSExpressionTree *expTree = NULL;
03575         // now initialize starts and variable index maps 
03576         for ( i = 0; i < iNumJacRowStarts; i++){                        
03577                 m_miJacStart [ i ] = 0;
03578                 // map the variables  in the nonlinear rows
03579                 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
03580                         // the following is equivalent to  m_treeRoot->getVariableIndexMap( i);
03581                         m_mapExpressionTreesMod[ i]->getVariableIndiciesMap();
03582                         
03583                 }
03584         }
03585         int loopLimit =  getConstraintNumber();
03586         // only execute the following code if there are linear constraint coefficients
03587         // determine the number of terms in constraint with constant partial derivative
03588         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
03589                 // i is indexing rows (constrains) and j is indexing column numbers 
03590                 for (i = 0; i < loopLimit; i++){
03591                         m_miJacNumConTerms[ i] = 0;
03592                         for (j = start[i]; j < start[ i + 1 ]; j++){
03593                                 // determine if variable index[j] appears in the Expression Tree for row i
03594                                 // if we pass if test below then variable i is in the expresssion tree and we add
03595                                 // the linear term to the expession tree
03596                                 if( (m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ) &&
03597                                         ( (*m_mapExpressionTreesMod[ i]->mapVarIdx).find( index[ j]) != (*m_mapExpressionTreesMod[ i]->mapVarIdx).end()) ){
03598                                         // variable index[ j] appears in the expression tree for row i
03599                                         // add the coefficient corresponding to variable index[j] in row i to the expression tree       
03600                                         // define a new OSnLVariable and OSnLnodePlus 
03601                                         if(value[ j] > 0 || value[j] < 0){
03602                                                 nlNodeVariable = new OSnLNodeVariable();
03603                                                 nlNodeVariable->coef = value[ j];
03604                                                 nlNodeVariable->idx = index[ j];
03605                                                 nlNodePlus = new OSnLNodePlus();
03606                                                 nlNodePlus->m_mChildren[ 0] = m_mapExpressionTreesMod[ i ]->m_treeRoot;
03607                                                 nlNodePlus->m_mChildren[ 1] = nlNodeVariable;
03608                                                 //expTree = new OSExpressionTree();
03609                                                 //expTree->m_treeRoot = nlNodePlus ;
03610                                                 //expTree->mapVarIdx = m_mapExpressionTreesMod[ i]->mapVarIdx;
03611                                                 //m_mapExpressionTreesMod[ i ]  = expTree;      
03612                                                 m_mapExpressionTreesMod[ i ]->m_treeRoot = nlNodePlus;
03613                                         }
03614                                 }
03615                                 else{ 
03616                                         //the partial derivative of variable j in row i is constant
03617                                         m_miJacNumConTerms[ i]++;
03618                                 }                               
03619                         }
03620                 }
03621         }
03622         //
03623         m_miJacStart[0] = 0;
03624         for (i = 1; i < iNumJacRowStarts; i++ ){
03625                 if( m_mapExpressionTreesMod.find( i - 1) != m_mapExpressionTreesMod.end() ){
03626                         m_miJacStart[i] = m_miJacStart[i - 1] + (m_miJacNumConTerms[ i - 1] + (*m_mapExpressionTreesMod[ i - 1]->mapVarIdx).size() );
03627                 }
03628                 else{
03629                         m_miJacStart[i] = m_miJacStart[i - 1] + m_miJacNumConTerms[ i - 1];
03630                 }       
03631         }
03632         // we know how many constant terms and size of arrays
03633         // dimension miIndex and mdValue here
03634         m_iJacValueSize =       m_miJacStart[ iNumJacRowStarts - 1];
03635         m_miJacIndex = new int[  m_iJacValueSize];
03636         m_mdJacValue = new double[ m_iJacValueSize ];
03637         // now loop again and put in values and indicies
03638         // first put in the constant terms
03639         if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
03640                 for (i = 0; i < loopLimit; i++){
03641                         k = 0;
03642                         for (j = start[i]; j < start[ i + 1 ]; j++){
03643                                 if( (m_mapExpressionTreesMod.find( i) == m_mapExpressionTreesMod.end() ) ||
03644                                         ( (*m_mapExpressionTreesMod[ i]->mapVarIdx).find( index[ j]) == (*m_mapExpressionTreesMod[ i]->mapVarIdx).end()) ){
03645                                                 m_miJacIndex[ m_miJacStart[i] + k ] = index[ j];
03646                                                 m_mdJacValue[ m_miJacStart[i] + k ] = value[ j];
03647                                                 k++;
03648                                         }                       
03649                         }
03650                 }
03651         }
03652         // put in terms from the modified nonlinear expression tree
03653         for (i = 0; i < loopLimit; i++ ){
03654                 k = m_miJacStart[i] + m_miJacNumConTerms[ i ];
03655                 // if the row is in the list of expression trees read in idicies and values
03656                 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
03657                         for(posVarIdx = (*m_mapExpressionTreesMod[ i]->mapVarIdx).begin(); posVarIdx 
03658                         != (*m_mapExpressionTreesMod[ i]->mapVarIdx).end(); ++posVarIdx){
03659                                 m_miJacIndex[ k] = posVarIdx->first;
03660                                 m_mdJacValue[ k] = 0;
03661                                 k++;
03662                         }
03663                 }
03664         }
03665         #ifdef DEBUG
03666         std::cout << "HERE ARE ROW STARTS:" << std::endl;
03667         for (i = 0; i < iNumJacRowStarts; i++ ){
03668                 std::cout <<  m_miJacStart[ i] << "  "; 
03669         }
03670         std::cout << std::endl << std::endl;
03671         std::cout << "HERE ARE VARIABLE INDICIES:" << std::endl;
03672         for (i = 0; i < m_miJacStart[ iNumJacRowStarts - 1]; i++ ){
03673                 std::cout <<  m_miJacIndex[ i] << "  "; 
03674         }
03675         std::cout << std::endl << std::endl;
03676         std::cout << "HERE ARE VALUES:" << std::endl;
03677         for (i = 0; i < m_miJacStart[ iNumJacRowStarts - 1]; i++ ){
03678                 std::cout <<  m_mdJacValue[ i] << "  "; 
03679         }
03680         std::cout << std::endl << std::endl;
03681 
03682         std::cout << "HERE ARE NUMBER OF CONSTANT TERMS:" << std::endl;
03683         for (i = 0; i < iNumJacRowStarts - 1; i++ ){
03684                 std::cout <<  m_miJacNumConTerms[ i ] << "  ";  
03685         }
03686         std::cout << std::endl << std::endl;
03687         #endif
03688         return true;
03689 }//getSparseJacobianFromRowMajor
03690 
03691 OSExpressionTree* OSInstance::getLagrangianExpTree( ){
03692         if( m_bLagrangianExpTreeCreated == true) return m_LagrangianExpTree;
03693         // we calculate the Lagrangian for all the objectives and constraints
03694         // with nonlinear terms
03695         // first initialize everything for nonlinear work
03696         if(m_bSparseJacobianCalculated == false) getJacobianSparsityPattern( ); 
03697         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03698         OSnLNodeTimes* nlNodeTimes = NULL;
03699         OSnLNodeVariable* nlNodeVariable = NULL;
03700         OSnLNodeSum* nlNodeSum = NULL;
03701         int numChildren = 0;
03702         int rowIdx;
03703         // create the sum node
03704         nlNodeSum = new OSnLNodeSum();
03705         nlNodeSum->inumberOfChildren = m_mapExpressionTreesMod.size();
03706         nlNodeSum->m_mChildren = new OSnLNode*[ nlNodeSum->inumberOfChildren];
03707         // create and expression tree for the sum node
03708         m_LagrangianExpTree = new OSExpressionTree();
03709         m_LagrangianExpTree->m_treeRoot = nlNodeSum;
03710         // now create the children of the sum node
03711         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03712                 // this variable is the Lagrange multiplier
03713                 nlNodeVariable = new OSnLNodeVariable();
03714                 nlNodeVariable->coef = 1.;
03715                 // get the correct index --
03716                 // for rowIdx = 0, ..., m - 1 set idx = numVar + rowIdx
03717                 rowIdx = posMapExpTree->first;
03718                 if(rowIdx >= 0){
03719                         nlNodeVariable->idx = instanceData->variables->numberOfVariables + rowIdx;
03720                 }
03721                 else{
03722                         nlNodeVariable->idx = instanceData->variables->numberOfVariables + 
03723                         instanceData->constraints->numberOfConstraints + (abs(rowIdx) - 1);
03724                 }
03725                 // now create a times multiply the new variable times the root of the expression tree
03726                 nlNodeTimes = new OSnLNodeTimes();
03727                 nlNodeTimes->m_mChildren[ 0] = nlNodeVariable;
03728                 nlNodeTimes->m_mChildren[ 1] = m_mapExpressionTreesMod[ posMapExpTree->first ]->m_treeRoot;     
03729                 // the times node is the new child
03730                 nlNodeSum->m_mChildren[ numChildren] = nlNodeTimes;
03731                 numChildren++;
03732         }       
03733         // get a variable index map for the expression tree
03734         m_LagrangianExpTree->getVariableIndiciesMap();
03735         // print out the XML for this puppy
03736         //std::cout << "Here comes the Lagrangian Tree" << std::endl;
03737         //std::cout << m_LagrangianExpTree->m_treeRoot->getNonlinearExpressionInXML() << std::endl;
03738         //
03739         m_bLagrangianExpTreeCreated = true;
03740         return m_LagrangianExpTree;
03741 }//getLagrangianExpTree
03742 
03743 std::map<int, int> OSInstance::getAllNonlinearVariablesIndexMap( ){
03744         if(m_bAllNonlinearVariablesIndex == true) return m_mapAllNonlinearVariablesIndex;
03745         //loop over the map of expression tree and get a unique listing of all variables
03746         // put these in the map m_mapAllNonlinearVariablesIndex
03747         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03748         std::map<int, int>::iterator posVarIdx;
03749         OSExpressionTree *expTree;
03750         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03751                 // get the index map for the expression tree
03752                 expTree = posMapExpTree->second;
03753                 if(expTree->m_bIndexMapGenerated == false)expTree->getVariableIndiciesMap();
03754                 for(posVarIdx = (*expTree->mapVarIdx).begin(); posVarIdx != (*expTree->mapVarIdx).end(); ++posVarIdx){
03755                         if( m_mapAllNonlinearVariablesIndex.find( posVarIdx->first) == m_mapAllNonlinearVariablesIndex.end() ){
03756                         // add the variable to the Lagragian map
03757                         m_mapAllNonlinearVariablesIndex[ posVarIdx->first] = 1;
03758                         }
03759                 }
03760         }
03761         m_miNonLinearVarsReverseMap = new int[m_mapAllNonlinearVariablesIndex.size()];
03762         // now order appropriately
03763         int kount = 0;
03764         //std::cout << "HERE IS THE LAGRANGIANN VARIABLE MAPPING" << std::endl;
03765         for(posVarIdx = m_mapAllNonlinearVariablesIndex.begin(); posVarIdx !=m_mapAllNonlinearVariablesIndex.end(); ++posVarIdx){
03766                 posVarIdx->second = kount;
03767                 m_miNonLinearVarsReverseMap[ kount] = posVarIdx->first;
03768                 kount++;
03769                 #ifdef DEBUG
03770                 std::cout <<  "POSITION FIRST =  "  << posVarIdx->first ;
03771                 std::cout <<  "   POSITION SECOND = "  << posVarIdx->second << std::endl;
03772                 #endif
03773         }
03774         m_iNumberOfNonlinearVariables = kount;
03775         //std::cout <<  "NUMBER OF NONLINEAR VARIABLES =  "  << kount ;
03776         m_bAllNonlinearVariablesIndex = true;
03777         return m_mapAllNonlinearVariablesIndex;
03778 }//getAllNonlinearVariablesIndexMap     
03779 
03780 SparseHessianMatrix* OSInstance::getLagrangianHessianSparsityPattern( ){
03781         // fill in the nonzeros in the sparse Hessian
03782         if( m_bLagrangianSparseHessianCreated == true) return m_LagrangianSparseHessian;
03783         if( m_iNumberOfNonlinearVariables == 0) return NULL;
03784         if( m_binitForAlgDiff == false ) initForAlgDiff();
03785         unsigned int i = 0;
03786         unsigned int j;
03787         int numNonz = 0;
03788         // Create the CppAD function if necessary
03789         //
03790         std::vector<double> vx;
03791         std::map<int, int>::iterator posMap1, posMap2;  
03792         if( (m_bOSADFunIsCreated == false || m_bCppADMustReTape == true )  && (m_mapExpressionTreesMod.size() > 0) ) {
03793                 for(posMap1 = m_mapAllNonlinearVariablesIndex.begin(); posMap1 != m_mapAllNonlinearVariablesIndex.end(); ++posMap1){
03794                         vx.push_back( 1.0) ;
03795                 }
03796                 createOSADFun( vx);
03797         }
03798         
03799         
03800         /*
03801         //
03802         // Use CppAD to do a forward sparsity calculation
03803         std::vector<bool> r(m_iNumberOfNonlinearVariables * m_iNumberOfNonlinearVariables);
03804         unsigned int j;
03805         for(i = 0; i < m_iNumberOfNonlinearVariables; i++) { 
03806                 for(j = 0; j < m_iNumberOfNonlinearVariables; j++)
03807                         r[ i * m_iNumberOfNonlinearVariables + j ] = false;
03808                         r[ i * m_iNumberOfNonlinearVariables + i] = true;
03809         }       
03810         // compute sparsity pattern for J(x) = F^{(1)} (x)
03811         (*Fad).ForSparseJac(m_iNumberOfNonlinearVariables, r);
03812         //
03813         //now the second derivative
03814         unsigned int m = m_mapExpressionTreesMod.size();
03815         std::vector<bool> e( m);
03816         //Vector s(m);
03817         for(i = 0; i < m; i++) e[i] = true;
03818         std::cout << "Computing Sparse Hessian" << std::endl;
03819         //m_vbLagHessNonz holds the sparsity pattern Lagrangian of the Hessian
03820         m_vbLagHessNonz = (*Fad).RevSparseHes(m_iNumberOfNonlinearVariables, e);
03821         for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
03822                 //std::cout << "Row " << i << "  of Hessian " << std::endl;
03823                 for(j = i; j < m_iNumberOfNonlinearVariables; j++){
03824                         if(m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j]  == true) numNonz++;
03825                         //std::cout << m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j] <<  "  " ;
03826                 }
03827                 //std::cout << std::endl;
03828         }
03829         
03830         */
03831         // note this function call also fills in m_vbLagHessNonz
03832         numNonz = getADSparsityHessian();
03833         
03834         //std::cout << "Lagrangian Hessian Nonzeros = " << numNonz << std::endl;
03835         i = 0;
03836         // now that we have the dimension create SparseHessianMatrix (upper triangular)
03837         m_LagrangianSparseHessian = new SparseHessianMatrix();
03838         m_LagrangianSparseHessian->bDeleteArrays = true;
03839         m_LagrangianSparseHessian->hessDimension = numNonz;
03840         //m_LagrangianSparseHessian->hessDimension = m_vbLagHessNonz.size();
03841         m_LagrangianSparseHessian->hessRowIdx = new int[m_LagrangianSparseHessian->hessDimension];
03842         m_LagrangianSparseHessian->hessColIdx = new int[m_LagrangianSparseHessian->hessDimension];
03843         m_LagrangianSparseHessian->hessValues = new double[m_LagrangianSparseHessian->hessDimension];
03844         //std::cout << "HESSIAN DIMENSION = " << m_LagrangianSparseHessian->hessDimension << std::endl;
03845         numNonz = 0;
03846         for(posMap1 = m_mapAllNonlinearVariablesIndex.begin(); posMap1 != m_mapAllNonlinearVariablesIndex.end(); ++posMap1){
03847                 //std::cout << "posMap1->first  " << posMap1->first << std::endl;
03848                 j = i;
03849                 for(posMap2 = posMap1; posMap2 != m_mapAllNonlinearVariablesIndex.end(); ++posMap2){
03850                         //std::cout << "posMap2->first  " << posMap2->first << std::endl;
03851                         if(m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j] == true){
03852                                 *(m_LagrangianSparseHessian->hessRowIdx + numNonz) = posMap1->first;
03853                                 *(m_LagrangianSparseHessian->hessColIdx + numNonz) = posMap2->first;
03854                                 numNonz++;
03855                         }
03856                         //std::cout << m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j] <<  "  " << std::endl;
03857                         j++;
03858                 }
03859                 i++;
03860         }
03861         #ifdef DEBUG
03862         std::cout << "HESSIAN SPARSITY PATTERN" << std::endl;
03863         int kj;
03864         for(kj = 0; kj < m_LagrangianSparseHessian->hessDimension; kj++){
03865                 std::cout <<  "Row Index = " << *(m_LagrangianSparseHessian->hessRowIdx + kj) << std::endl;
03866                 std::cout <<  "Column Index = " << *(m_LagrangianSparseHessian->hessColIdx + kj) << std::endl;
03867         }
03868         #endif
03869         //
03870         m_bLagrangianSparseHessianCreated = true;
03871         return m_LagrangianSparseHessian;
03872 }//getLagrangianHessianSparsityPattern
03873 
03874 
03875 void OSInstance::duplicateExpressionTreesMap(){
03876         //std::cout << "I AM IN DUPLICATE EXPRSSION TREES MAP" << std::endl;
03877         // we do this so that we can keep the integrity of the original formulation
03878         if(m_bDuplicateExpressionTreesMap == false){ 
03879                 // first make sure the map was created
03880                 if( m_bProcessExpressionTrees == false) getAllNonlinearExpressionTrees();
03881                 m_mapExpressionTreesMod = m_mapExpressionTrees;
03882                 m_bDuplicateExpressionTreesMap = true;
03883                 return;
03884         }
03885         else{
03886                 return;
03887         }
03888 }//duplicateExpressionTreesMap
03889 
03890 
03891 
03892 bool OSInstance::getIterateResults( double *x, double *objLambda, double* conMultipliers, 
03893                 bool new_x, int highestOrder){
03894         try{ 
03895                 if( m_binitForAlgDiff == false) initForAlgDiff();
03896                 std::map<int, int>::iterator posVarIndexMap;
03897                 
03898                 if(new_x == true){
03899                         if( m_vdX.size() > 0) m_vdX.clear();
03900                         for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
03901                                 m_vdX.push_back( x[ posVarIndexMap->first]) ;
03902                         }
03903                         if( (m_bOSADFunIsCreated == false || m_bCppADMustReTape == true )  && (m_mapExpressionTreesMod.size() > 0) ) {
03904                                 createOSADFun( m_vdX);
03905                         }       
03906                 }
03907                 else{ // make sure vector not empty -- this could happen if we have linear obj and nonlinear constraints
03908                         if( m_vdX.size() == 0) {
03909                                 for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
03910                                         m_vdX.push_back( x[ posVarIndexMap->first]) ;
03911                                 }
03912                                 if( (m_bOSADFunIsCreated == false || m_bCppADMustReTape == true )  && (m_mapExpressionTreesMod.size() > 0) ) {
03913                                         createOSADFun( m_vdX);
03914                                 }       
03915                         }
03916                 }
03917                 switch( highestOrder){          
03918                         case 0: 
03919                                 if(new_x == true || m_iHighestOrderEvaluated < 0){      
03920                                         if(bUseExpTreeForFunEval == true){
03921                                                 calculateAllConstraintFunctionValues( x, new_x);
03922                                                 calculateAllObjectiveFunctionValues( x, new_x);
03923                                         }
03924                                         else{
03925                                                 getZeroOrderResults(x, objLambda, conMultipliers);
03926                                         }
03927 
03928                                 }
03929                                 break;  
03930                         case 1:
03931                                 if(new_x == true || m_iHighestOrderEvaluated < 0)       
03932                                         getZeroOrderResults(x, objLambda, conMultipliers);
03933                                 if(new_x == true || m_iHighestOrderEvaluated < 1)       
03934                                         getFirstOrderResults(x, objLambda, conMultipliers);
03935                                 break;
03936                         case 2: 
03937                                 if(new_x == true || m_iHighestOrderEvaluated < 0)       
03938                                         getZeroOrderResults(x, objLambda, conMultipliers);
03939                                 if(new_x == true || m_iHighestOrderEvaluated < 2)       
03940                                         getSecondOrderResults(x, objLambda, conMultipliers);
03941                                 break;
03942                         default:
03943                                 throw ErrorClass("Derivative should be order 0, 1, or 2");      
03944                 }//end switch
03945                 return true;
03946         }
03947         catch(const ErrorClass& eclass){
03948                 throw ErrorClass( eclass.errormsg);
03949         }  
03950 }//end getIterateResults
03951 
03952 
03953 bool OSInstance::getZeroOrderResults(double *x, double *objLambda, double *conMultipliers){
03954         try{ 
03955                 // initialize everything
03956                 int i, j, rowNum, objNum;
03957                 if( m_mapExpressionTreesMod.size() > 0){
03958                         m_vdYval = this->forwardAD(0, m_vdX);   
03959                 }
03960                 // now get all function and constraint values using forward result
03961                 for(rowNum = 0; rowNum < m_iConstraintNumber; rowNum++){
03962                         m_mdConstraintFunctionValues[ rowNum] = 0.0;
03963                         if( m_mapExpressionTreesMod.find( rowNum) != m_mapExpressionTreesMod.end() ){
03964                                 m_mdConstraintFunctionValues[ rowNum] = m_vdYval[  m_mapOSADFunRangeIndex[ rowNum]];
03965                         }
03966                         // now the linear part
03967                         // be careful, loop over only the constant terms in sparseJacMatrix
03968                         i = m_sparseJacMatrix->starts[ rowNum];
03969                         j = m_sparseJacMatrix->starts[ rowNum + 1 ];
03970                         while ( (i - m_sparseJacMatrix->starts[ rowNum])  < m_sparseJacMatrix->conVals[ rowNum] ){
03971                                 m_mdConstraintFunctionValues[ rowNum] += m_sparseJacMatrix->values[ i]*x[ m_sparseJacMatrix->indexes[ i] ];
03972                                 i++;
03973                         }       
03974                         // add in the constraint function constant
03975                         m_mdConstraintFunctionValues[ rowNum] += m_mdConstraintConstants[ rowNum ];
03976                         #ifdef DEBUG
03977                         std::cout << "Constraint " <<  rowNum << " function value =  " << m_mdConstraintFunctionValues[ rowNum ] << std::endl;
03978                         #endif
03979                 }
03980                 // now get the objective function values from the forward result
03981                 for(objNum = 0; objNum < m_iObjectiveNumber; objNum++){
03982                         m_mdObjectiveFunctionValues[ objNum] = 0.0;
03983                         if( m_mapExpressionTreesMod.find( -objNum -1) != m_mapExpressionTreesMod.end() ){
03984                                 m_mdObjectiveFunctionValues[ objNum] = m_vdYval[ objNum];
03985                         }
03986                         for(i = 0; i < m_iVariableNumber; i++){
03987                                 m_mdObjectiveFunctionValues[ objNum] += m_mmdDenseObjectiveCoefficients[ objNum][i]*x[ i];
03988                         }
03989                         #ifdef DEBUG
03990                         std::cout << "Objective " << objNum << " function value =  " << m_mdObjectiveFunctionValues[ objNum] << std::endl;
03991                         #endif
03992                 }
03993         return true;
03994         }//end try
03995         catch(const ErrorClass& eclass){
03996                 throw ErrorClass( eclass.errormsg);
03997         }  
03998 }//end getZeroOrderResults
03999 
04000 
04001 
04002 bool OSInstance::getFirstOrderResults(double *x, double *objLambda, double *conMultipliers){
04003         try{
04004                 // initialize everything
04005                 unsigned int i, j;
04006                 int rowNum,  jacIndex;
04007                 unsigned int jstart, jend;
04008                 int idx;
04009                 OSExpressionTree *expTree = NULL;
04010                 int domainIdx = 0;      
04011                 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
04012                 std::map<int, int>::iterator posVarIdx;
04013                 
04019                 if(m_iNumberOfNonlinearVariables >= m_mapExpressionTreesMod.size() ){
04020                         // calculate the gradient by doing a reverse sweep over each row
04021                         // loop over the constraints that have a nonlinear term and get their gradients
04022                         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
04023                                 idx = posMapExpTree->first;
04024                                 // we are considering only constraints, not objective function
04025                                 if(idx >= 0){
04026                                         m_vdRangeUnitVec[ domainIdx] = 1.;
04027                                         m_mapExpressionTreesMod[ idx]->getVariableIndiciesMap(); 
04028                                         m_vdYjacval = this->reverseAD(1, m_vdRangeUnitVec);
04029                                         // check size
04030                                         jstart = m_miJacStart[ idx] + m_miJacNumConTerms[ idx];
04031                                         jend = m_miJacStart[ idx + 1 ];
04032                                         if( (*m_mapExpressionTreesMod[ idx]->mapVarIdx).size() != (jend - jstart)) throw 
04033                                         ErrorClass("number of partials not consistent");
04034                                         j = 0;
04035                                         jacIndex = 0;
04036                                         for(posVarIdx = m_mapAllNonlinearVariablesIndex.begin(); posVarIdx 
04037                                                 != m_mapAllNonlinearVariablesIndex.end(); ++posVarIdx){
04038                                                 //std::cout << "Constraint Function Jacobian Values" << "For Constraint  " << idx  << std::endl;
04039                                                 //std::cout << "Jac Val for index " << posVarIdx->first  << " = " << m_vdYjacval[ jacIndex] << std::endl;
04040                                                 //if(m_miJacIndex[ jstart] != posVarIdx->first) throw ErrorClass("error calculating Jacobian matrix");
04041                                                 // we are working with variable posVarIdx->first in the original variable space
04042                                                 // we need to see which variable this is in the individual constraint map
04043                                                 if( (*m_mapExpressionTreesMod[ idx]->mapVarIdx).find( posVarIdx->first) != (*m_mapExpressionTreesMod[ idx]->mapVarIdx).end()){
04044                                                         m_mdJacValue[ jstart] = m_vdYjacval[ jacIndex];
04045                                                         jstart++;
04046                                                         j++;
04047                                                 }
04048                                                 jacIndex++;
04049                                         }
04050                                         
04051                                         m_vdRangeUnitVec[ domainIdx] = 0.;
04052                                         domainIdx++;
04053                                 }
04054                                 else{    // we have an objective function
04055                                         
04056 
04057                                         domainIdx++;
04058                                 }
04059                         }
04060                 }
04061                 else{  
04062                         // calculate the gradients using a forward sweep over all the variables.                
04063                         for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
04064                                 m_vdDomainUnitVec[i] = 1.;     
04065                                 rowNum = 0;
04066                                 if( m_mapExpressionTreesMod.size() > 0){          
04067                                         m_vdYjacval = this->forwardAD(1, m_vdDomainUnitVec);
04068                                 } 
04069                                 // fill in Jacobian here, we have column i 
04070                                 // start Jacobian calculation
04071                                 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
04072                                         idx = posMapExpTree->first;
04073                                         // we are considering only constraints, not objective function
04074                                         if(idx >= 0){
04075                                                 //figure out original variable this corresponds to
04076                                                 //then use (*m_mapExpressionTreesMod[ idx]->mapVarIdx) to figure out which variable it is within row idx
04077                                                 //m_mapAllNonlinearVariablesIndex
04078                                                 //std::cout << "This is the following variable in the expression tree  " <<  (*m_mapExpressionTreesMod[ idx]->mapVarIdx)[ m_miNonLinearVarsReverseMap[ i]]<< std::endl;                         
04079                                                 expTree = m_mapExpressionTreesMod[ idx];                
04080                                                 if( (*expTree->mapVarIdx).find( m_miNonLinearVarsReverseMap[ i]) != (*expTree->mapVarIdx).end()  ){             
04081                                                         jacIndex = (*m_mapExpressionTreesMod[ idx]->mapVarIdx)[ m_miNonLinearVarsReverseMap[ i]];
04082                                                         jstart = m_miJacStart[ idx] + m_miJacNumConTerms[ idx];
04083                                                         // kipp change 1 to number of objective functions
04084                                                         m_mdJacValue[ jstart + jacIndex] = m_vdYjacval[m_iObjectiveNumberNonlinear + rowNum];
04085                                                 }
04086                                                 rowNum++;
04087                                         }//end Jacobian calculation
04088                         }                       
04089                         //
04090                         m_vdDomainUnitVec[i] = 0.;
04091                         }
04092                 }
04093                 #ifdef DEBUG
04094                 int k;
04095                 std::cout  << "JACOBIAN DATA " << std::endl;
04096                 for(idx = 0; idx < m_iConstraintNumber; idx++){
04097                         for(k = *(m_sparseJacMatrix->starts + idx); k < *(m_sparseJacMatrix->starts + idx + 1); k++){
04098                                 std::cout << "row idx = " << idx <<  "  col idx = "<< *(m_sparseJacMatrix->indexes + k)
04099                                 << " value = " << *(m_sparseJacMatrix->values + k) << std::endl;
04100                         }
04101                 }
04102                 #endif
04103                 return true;
04104         }//end try
04105         catch(const ErrorClass& eclass){
04106                 throw ErrorClass( eclass.errormsg);
04107         } 
04108 }// end getFirstOrderResults
04109                         
04110 
04111 bool OSInstance::getSecondOrderResults(double *x, double *objLambda, double *conMultipliers){
04112         try{
04113                 // initialize everything
04114                 unsigned int i, j;
04115                 int rowNum,  jacIndex;
04116                 int jstart,  idx;
04117                 OSExpressionTree *expTree = NULL;
04118                 int hessValuesIdx = 0;  
04119                 if( m_bLagrangianSparseHessianCreated == false) getLagrangianHessianSparsityPattern( );
04120                 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
04121                 std::map<int, int>::iterator posVarIndexMap;
04122                 if( objLambda == NULL) throw ErrorClass("must have a multiplier for the objective function even if zero when calling getSecondOrderResults");
04123                 if( conMultipliers == NULL) throw ErrorClass("cannot have a null vector of lagrange multipliers when calling getSecondOrderResults -- okay if  zero");
04124                 if( m_vdLambda.size() > 0) m_vdLambda.clear();
04125                 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){  
04126                         if( posMapExpTree->first >= 0){
04127                                 m_vdLambda.push_back( conMultipliers[ posMapExpTree->first]);
04128                         }
04129                         else{
04130                                 // kipp correct when there is more than one obj
04131                                 m_vdLambda.push_back( objLambda[ abs(posMapExpTree->first) - 1] );
04132                         }
04133                 }
04134                 for(i = 0; i < m_iNumberOfNonlinearVariables; i++){ 
04135                         m_vdDomainUnitVec[i] = 1.;     
04136                         rowNum = 0;
04137                         if( m_mapExpressionTreesMod.size() > 0){          
04138                                 m_vdYjacval = this->forwardAD(1, m_vdDomainUnitVec);
04139                         } 
04140                         // fill in Jacobian here, we have column i 
04141                         // start Jacobian calculation
04142                         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
04143                                 idx = posMapExpTree->first;
04144                                 // we are considering only constraints, not objective function
04145                                 if(idx >= 0){
04146                                         //figure out original variable this corresponds to
04147                                         //then use (*m_mapExpressionTreesMod[ idx]->mapVarIdx) to figure out which variable it is within row idx
04148                                         //m_mapAllNonlinearVariablesIndex
04149                                         //std::cout << "This is the following variable in the expression tree  " <<  (*m_mapExpressionTreesMod[ idx]->mapVarIdx)[ m_miNonLinearVarsReverseMap[ i]]<< std::endl;                         
04150                                         expTree = m_mapExpressionTreesMod[ idx];                
04151                                         if( (*expTree->mapVarIdx).find( m_miNonLinearVarsReverseMap[ i]) != (*expTree->mapVarIdx).end()  ){             
04152                                                 jacIndex = (*m_mapExpressionTreesMod[ idx]->mapVarIdx)[ m_miNonLinearVarsReverseMap[ i]];
04153                                                 jstart = m_miJacStart[ idx] + m_miJacNumConTerms[ idx];
04154                                                 m_mdJacValue[ jstart + jacIndex] = m_vdYjacval[m_iObjectiveNumberNonlinear + rowNum];
04155                                         }
04156                                         rowNum++;
04157                                 }//end Jacobian calculation
04158                                 else{
04159                                         // see if we have the objective function of interest
04160                                         //kipp fix if more than one obj
04161                                                 m_mmdObjGradient[  (abs( idx) - 1)][ m_miNonLinearVarsReverseMap[ i]] = m_vdYjacval[ (abs( idx) - 1)] + 
04162                                                 m_mmdDenseObjectiveCoefficients[  (abs( idx) - 1)][ m_miNonLinearVarsReverseMap[ i]];                                   
04163                         }//end Obj gradient calculation 
04164                 }                       
04165                 // now calculate the Hessian
04166                 if( m_mapExpressionTreesMod.size() > 0){   
04167                         m_vdw = reverseAD(2, m_vdLambda);   // derivative of partial
04168                 }
04169                 for(j = i; j < m_iNumberOfNonlinearVariables; j++){
04170                         if( m_vbLagHessNonz[i*m_iNumberOfNonlinearVariables + j] == true){
04171                                 m_LagrangianSparseHessian->hessValues[ hessValuesIdx] =  m_vdw[  j*2 + 1];
04172                                 #ifdef DEBUG
04173                                 std::cout << "reverse 2 " << m_LagrangianSparseHessian->hessValues[ hessValuesIdx] << std::endl;
04174                                 #endif
04175                                 hessValuesIdx++;
04176                         }
04177                 }
04178                 //
04179                 //
04180                 m_vdDomainUnitVec[i] = 0.;
04181         }
04182         #ifdef DEBUG
04183         int k;
04184         std::cout  << "JACOBIAN DATA " << std::endl;
04185         for(idx = 0; idx < m_iConstraintNumber; idx++){
04186                 for(k = *(m_sparseJacMatrix->starts + idx); k < *(m_sparseJacMatrix->starts + idx + 1); k++){
04187                         std::cout << "row idx = " << idx <<  "  col idx = "<< *(m_sparseJacMatrix->indexes + k)
04188                         << " value = " << *(m_sparseJacMatrix->values + k) << std::endl;
04189                 }
04190         }
04191         #endif
04192         return true;
04193         }//end try
04194         catch(const ErrorClass& eclass){
04195                 throw ErrorClass( eclass.errormsg);
04196         } 
04197 }// end getSecondOrderResults
04198 
04199 bool OSInstance::initForAlgDiff(){
04200         if( m_binitForAlgDiff == true ) return true;
04201         initializeNonLinearStructures( );
04202         initObjGradients();
04203         getAllNonlinearVariablesIndexMap( );
04204         //if(m_bSparseJacobianCalculated  == false) getJacobianSparsityPattern();
04205         //see if we need to retape 
04206         //loop over expression tree and see if one requires it
04207         std::map<int, OSExpressionTree*>::iterator posMapExpTree;
04208         for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
04209                 if(posMapExpTree->second->bADMustReTape == true) m_bCppADMustReTape = true;
04210         }                               
04211 
04212         #ifdef DEBUG
04213         std::cout << "RETAPE ==  " << m_bCppADMustReTape << std::endl;
04214         #endif
04215         unsigned int i;
04216         for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
04217                 m_vdDomainUnitVec.push_back( 0.0 );
04218         }
04219         for(i = 0; i < m_mapExpressionTreesMod.size(); i++){
04220                 m_vdRangeUnitVec.push_back( 0.0 );
04221         }
04222         m_binitForAlgDiff = true;
04223         //m_bSparseJacobianCalculated = true;
04224         //m_bProcessExpressionTrees = true;
04225         return true;
04226 }//end initForAlgDiff
04227 
04228 bool OSInstance::initObjGradients(){
04229         int i, j;
04230         int m, n;
04231         m = getObjectiveNumber();
04232         n = getVariableNumber();
04233         getDenseObjectiveCoefficients();
04234         m_mmdObjGradient = new double*[m];
04235         for(i = 0; i < m; i++){
04236                 m_mmdObjGradient[i] = new double[n];
04237                 for(j = 0; j < n; j++){
04238                         m_mmdObjGradient[i][j] =  m_mmdDenseObjectiveCoefficients[ i][j];
04239                         #ifdef DEBUG
04240                         std::cout << "m_mmdObjGradient[i][j] = " << m_mmdObjGradient[i][j]  << std::endl;
04241                         #endif
04242                 }
04243         }
04244         return true;
04245 }//end initObjGradients
04254 bool OSInstance::setTimeDomain(std::string format)
04255 {   if ((format != "stages") && (format != "interval") && (format != "none")) 
04256                 return false;
04257         if (instanceData->timeDomain == NULL)
04258         {       instanceData->timeDomain = new TimeDomain();
04259         }
04260         if (format == "stages")
04261         {       if (instanceData->timeDomain->interval != NULL)
04262                 {       delete instanceData->timeDomain->interval;
04263                         instanceData->timeDomain->interval = NULL;
04264                 }
04265                 if (instanceData->timeDomain->stages == NULL)
04266                         instanceData->timeDomain->stages = new TimeDomainStages();
04267                 m_sTimeDomainFormat = format;
04268         }
04269         if (format == "interval")
04270         {       if (instanceData->timeDomain->stages != NULL)
04271                 {       delete instanceData->timeDomain->stages;
04272                         instanceData->timeDomain->stages = NULL;
04273                 }
04274                 if (instanceData->timeDomain->interval == NULL)
04275                         instanceData->timeDomain->interval = new TimeDomainInterval();
04276                 m_sTimeDomainFormat = format;
04277         }
04278         if (format == "none")
04279         {       if (instanceData->timeDomain->stages != NULL)
04280                 {       delete instanceData->timeDomain->stages;
04281                         instanceData->timeDomain->stages = NULL;
04282                 }
04283                 if (instanceData->timeDomain->interval != NULL)
04284                 {       delete instanceData->timeDomain->interval;
04285                         instanceData->timeDomain->interval = NULL;
04286                 }
04287                 m_sTimeDomainFormat = "";
04288         }
04289         return true;
04290 } //end setTimeDomain
04291 
04295 bool OSInstance::setTimeDomainStages(int number, std::string *names)
04296 {   if (instanceData->timeDomain == NULL)
04297                 instanceData->timeDomain = new TimeDomain();
04298         if (instanceData->timeDomain->interval != NULL)
04299                 return false;
04300         if (instanceData->timeDomain->stages == NULL)
04301         {       instanceData->timeDomain->stages = new TimeDomainStages;
04302         }
04303         else
04304         {       if (instanceData->timeDomain->stages->numberOfStages != number)
04305                 {       for (int i = 0; i < instanceData->timeDomain->stages->numberOfStages; i++)
04306                         {       if (instanceData->timeDomain->stages->stage[i]->variables != NULL)
04307                                 {       delete instanceData->timeDomain->stages->stage[i]->variables;
04308                                         instanceData->timeDomain->stages->stage[i]->variables = NULL;
04309                                 }
04310                                 if (instanceData->timeDomain->stages->stage[i]->constraints != NULL)
04311                                 {       delete instanceData->timeDomain->stages->stage[i]->constraints;
04312                                         instanceData->timeDomain->stages->stage[i]->constraints = NULL;
04313                                 }
04314                                 if (instanceData->timeDomain->stages->stage[i]->objectives != NULL)
04315                                 {       delete instanceData->timeDomain->stages->stage[i]->objectives;
04316                                         instanceData->timeDomain->stages->stage[i]->objectives = NULL;
04317                                 }
04318                                 delete instanceData->timeDomain->stages->stage[i];
04319                                 instanceData->timeDomain->stages->stage[i] = NULL;
04320                         }
04321                         delete []instanceData->timeDomain->stages->stage;
04322                         instanceData->timeDomain->stages->stage = NULL;
04323                 }
04324         }
04325         if (number != 0 )
04326         {       if (instanceData->timeDomain->stages->stage == NULL)
04327                         instanceData->timeDomain->stages->stage = new TimeDomainStage*[number];
04328                 for (int i = 0; i < number; i++)
04329                 {       instanceData->timeDomain->stages->stage[i] = new TimeDomainStage();
04330                 }
04331                 instanceData->timeDomain->stages->numberOfStages = number;
04332         }
04333         for (int i = 0; i < number; i++)
04334         //initial or empty vars, cons, objectives and set default to all objectives
04335         {       if (instanceData->timeDomain->stages->stage[i]->variables != NULL)
04336                 {       delete instanceData->timeDomain->stages->stage[i]->variables;
04337                         instanceData->timeDomain->stages->stage[i]->variables = NULL;
04338                 }
04339                 instanceData->timeDomain->stages->stage[i]->variables = new TimeDomainStageVariables();
04340                 if (instanceData->timeDomain->stages->stage[i]->constraints != NULL)
04341                 {       delete instanceData->timeDomain->stages->stage[i]->constraints;
04342                         instanceData->timeDomain->stages->stage[i]->constraints = NULL;
04343                 }
04344                 instanceData->timeDomain->stages->stage[i]->constraints = new TimeDomainStageConstraints();
04345                 if (instanceData->timeDomain->stages->stage[i]->objectives != NULL)
04346                 {       delete instanceData->timeDomain->stages->stage[i]->objectives;
04347                         instanceData->timeDomain->stages->stage[i]->objectives = NULL;
04348                 }
04349                 instanceData->timeDomain->stages->stage[i]->objectives = new TimeDomainStageObjectives();
04350                 instanceData->timeDomain->stages->stage[i]->objectives->numberOfObjectives = instanceData->objectives->numberOfObjectives;
04351                 instanceData->timeDomain->stages->stage[i]->objectives->obj = new TimeDomainStageObj*[instanceData->objectives->numberOfObjectives];
04352                 for (int j = 0; j < instanceData->objectives->numberOfObjectives; j++)
04353                 {       instanceData->timeDomain->stages->stage[i]->objectives->obj[j] = new TimeDomainStageObj();
04354                         instanceData->timeDomain->stages->stage[i]->objectives->obj[j]->idx = -(j+1);
04355                 }
04356                 if (names != NULL)
04357                         instanceData->timeDomain->stages->stage[i]->name = names[i];
04358         }
04359         return true;
04360 } //end setTimeDomainStages
04361 
04365 bool OSInstance::setTimeDomainStageVariablesOrdered(int numberOfStages, int *numberOfVariables, int *startIdx)
04366 {       if (instanceData->timeDomain == NULL)
04367                 instanceData->timeDomain = new TimeDomain();
04368         if (instanceData->timeDomain->interval != NULL)
04369                 return false;
04370         if (instanceData->timeDomain->stages == NULL)
04371                 instanceData->timeDomain->stages = new TimeDomainStages();      
04372         if (instanceData->timeDomain->stages != NULL)
04373         {       if ((instanceData->timeDomain->stages->numberOfStages != numberOfStages) &&
04374                         (instanceData->timeDomain->stages->numberOfStages != 0))
04375                          return false;
04376         }
04377         if (instanceData->timeDomain->stages->numberOfStages == 0)
04378         {       instanceData->timeDomain->stages->numberOfStages = numberOfStages;
04379                 if (instanceData->timeDomain->stages->stage == NULL)
04380                         instanceData->timeDomain->stages->stage = new TimeDomainStage*[numberOfStages];
04381                 for (int i = 0; i < numberOfStages; i++)
04382                         instanceData->timeDomain->stages->stage[i] = new TimeDomainStage();
04383         }
04384         int checksum = 0;
04385         for (int i = 0; i < numberOfStages; i++)
04386         //initial or empty vars, cons, objectives and set default to all objectives
04387         {       if (instanceData->timeDomain->stages->stage[i]->variables != NULL)
04388                 {       delete instanceData->timeDomain->stages->stage[i]->variables;
04389                         instanceData->timeDomain->stages->stage[i]->variables = NULL;
04390                 }
04391                 instanceData->timeDomain->stages->stage[i]->variables = new TimeDomainStageVariables();
04392                 instanceData->timeDomain->stages->stage[i]->variables->startIdx = startIdx[i];
04393                 instanceData->timeDomain->stages->stage[i]->variables->numberOfVariables = numberOfVariables[i];
04394                 checksum += numberOfVariables[i];
04395         }
04396         return (checksum == instanceData->variables->numberOfVariables);
04397 } //end setTimeDomainVariablesOrdered
04398 
04402 bool OSInstance::setTimeDomainStageVariablesUnordered(int numberOfStages, int *numberOfVariables, int **varIndex)
04403 {   if (instanceData->timeDomain == NULL)
04404                 instanceData->timeDomain = new TimeDomain();
04405         if (instanceData->timeDomain->interval != NULL)
04406                 return false;
04407         if (instanceData->timeDomain->stages == NULL)
04408                 instanceData->timeDomain->stages = new TimeDomainStages();      
04409         if (instanceData->timeDomain->stages != NULL)
04410         {       if ((instanceData->timeDomain->stages->numberOfStages != numberOfStages) &&
04411                         (instanceData->timeDomain->stages->numberOfStages != 0))
04412                          return false;
04413         }
04414         if (instanceData->timeDomain->stages->numberOfStages == 0 )
04415         {       instanceData->timeDomain->stages->numberOfStages = numberOfStages;
04416                 if (instanceData->timeDomain->stages->stage == NULL)
04417                         instanceData->timeDomain->stages->stage = new TimeDomainStage*[numberOfStages];
04418                 for (int i = 0; i < numberOfStages; i++)
04419                         instanceData->timeDomain->stages->stage[i] = new TimeDomainStage();
04420         }
04421         int checksum = 0;
04422         for (int i = 0; i < numberOfStages; i++)
04423         //initial or empty vars, cons, objectives and set default to all objectives
04424         {       if (instanceData->timeDomain->stages->stage[i]->variables != NULL)
04425                 {       delete instanceData->timeDomain->stages->stage[i]->variables;
04426                         instanceData->timeDomain->stages->stage[i]->variables = NULL;
04427                 }
04428                 instanceData->timeDomain->stages->stage[i]->variables = new TimeDomainStageVariables();
04429                 instanceData->timeDomain->stages->stage[i]->variables->numberOfVariables = numberOfVariables[i];
04430                 instanceData->timeDomain->stages->stage[i]->variables->var = new TimeDomainStageVar*[numberOfVariables[i]];
04431                 for (int j = 0; j < numberOfVariables[i]; j++)
04432                 {       instanceData->timeDomain->stages->stage[i]->variables->var[j] = new TimeDomainStageVar();
04433                         instanceData->timeDomain->stages->stage[i]->variables->var[j]->idx = varIndex[i][j];
04434                 }
04435                 checksum += numberOfVariables[i];
04436         }
04437         if (checksum != instanceData->variables->numberOfVariables) return false;
04438         int *checkvar = new int[instanceData->variables->numberOfVariables];
04439         for (int j = 0; j < instanceData->variables->numberOfVariables; j++)
04440                 checkvar[j] = -1;
04441         int k;
04442         for (int i = 0; i < numberOfStages; i++)
04443                 for (int j = 0; j < instanceData->timeDomain->stages->stage[i]->variables->numberOfVariables; j++)
04444                 {       k = instanceData->timeDomain->stages->stage[i]->variables->var[j]->idx;
04445                         if (checkvar[k] != -1) 
04446                         {       delete [] checkvar;
04447                                 checkvar = NULL;
04448                                 return false;
04449                         }
04450                         checkvar[k] = instanceData->timeDomain->stages->stage[i]->variables->var[j]->idx;
04451                 }
04452         delete [] checkvar;
04453         checkvar = NULL;                        
04454         return true;
04455 } //end setTimeDomainVariablesUnordered
04456 
04460 bool OSInstance::setTimeDomainStageConstraintsOrdered(int numberOfStages, int *numberOfConstraints, int *startIdx)
04461 {   if (instanceData->timeDomain == NULL)
04462                 instanceData->timeDomain = new TimeDomain();
04463         if (instanceData->timeDomain->interval != NULL)
04464                 return false;
04465         if (instanceData->timeDomain->stages == NULL)
04466                 instanceData->timeDomain->stages = new TimeDomainStages();      
04467         if (instanceData->timeDomain->stages != NULL)
04468         {       if ((instanceData->timeDomain->stages->numberOfStages != numberOfStages) &&
04469                         (instanceData->timeDomain->stages->numberOfStages != 0))
04470                          return false;
04471         }
04472         if (instanceData->timeDomain->stages->numberOfStages == 0 )
04473         {       instanceData->timeDomain->stages->numberOfStages = numberOfStages;
04474                 if (instanceData->timeDomain->stages->stage == NULL)
04475                         instanceData->timeDomain->stages->stage = new TimeDomainStage*[numberOfStages];
04476                 for (int i = 0; i < numberOfStages; i++)
04477                         instanceData->timeDomain->stages->stage[i] = new TimeDomainStage();
04478         }
04479         int checksum = 0;
04480         for (int i = 0; i < numberOfStages; i++)
04481         //initial or empty vars, cons, objectives and set default to all objectives
04482         {       if (instanceData->timeDomain->stages->stage[i]->constraints != NULL)
04483                 {       delete instanceData->timeDomain->stages->stage[i]->constraints;
04484                         instanceData->timeDomain->stages->stage[i]->constraints = NULL;
04485                 }
04486                 instanceData->timeDomain->stages->stage[i]->constraints = new TimeDomainStageConstraints();
04487                 instanceData->timeDomain->stages->stage[i]->constraints->startIdx = startIdx[i];
04488                 instanceData->timeDomain->stages->stage[i]->constraints->numberOfConstraints = numberOfConstraints[i];
04489                 checksum += numberOfConstraints[i];
04490         }
04491         return (checksum == instanceData->constraints->numberOfConstraints);
04492 } // end of setTimeStageConstraintsOrdered
04493 
04497 bool OSInstance::setTimeDomainStageConstraintsUnordered(int numberOfStages, int *numberOfConstraints, int **conIndex)
04498 {   if (instanceData->timeDomain == NULL)
04499                 instanceData->timeDomain = new TimeDomain();
04500         if (instanceData->timeDomain->interval != NULL)
04501                 return false;
04502         if (instanceData->timeDomain->stages == NULL)
04503                 instanceData->timeDomain->stages = new TimeDomainStages();      
04504         if (instanceData->timeDomain->stages != NULL)
04505         {       if ((instanceData->timeDomain->stages->numberOfStages != numberOfStages) &&
04506                         (instanceData->timeDomain->stages->numberOfStages != 0))
04507                          return false;
04508         }
04509         if ( instanceData->timeDomain->stages->numberOfStages == 0 )
04510         {       instanceData->timeDomain->stages->numberOfStages = numberOfStages;
04511                 if (instanceData->timeDomain->stages->stage == NULL)
04512                         instanceData->timeDomain->stages->stage = new TimeDomainStage*[numberOfStages];
04513                 for (int i = 0; i < numberOfStages; i++)
04514                         instanceData->timeDomain->stages->stage[i] = new TimeDomainStage();
04515         }
04516         int checksum = 0;
04517         for (int i = 0; i < numberOfStages; i++)
04518         //initial or empty vars, cons, objectives and set default to all objectives
04519         {       if (instanceData->timeDomain->stages->stage[i]->constraints != NULL)
04520                 {       delete instanceData->timeDomain->stages->stage[i]->constraints;
04521                         instanceData->timeDomain->stages->stage[i]->constraints = NULL;
04522                 }
04523                 instanceData->timeDomain->stages->stage[i]->constraints = new TimeDomainStageConstraints();
04524                 instanceData->timeDomain->stages->stage[i]->constraints->numberOfConstraints = numberOfConstraints[i];
04525                 instanceData->timeDomain->stages->stage[i]->constraints->con = new TimeDomainStageCon*[numberOfConstraints[i]];
04526                 for (int j = 0; j < numberOfConstraints[i]; j++)
04527                 {       instanceData->timeDomain->stages->stage[i]->constraints->con[j] = new TimeDomainStageCon();
04528                         instanceData->timeDomain->stages->stage[i]->constraints->con[j]->idx = conIndex[i][j];
04529                 }
04530                 checksum += numberOfConstraints[i];
04531         }
04532         if (checksum != instanceData->constraints->numberOfConstraints) return false;
04533         int *checkvar = new int[instanceData->constraints->numberOfConstraints];
04534         for (int j = 0; j < instanceData->constraints->numberOfConstraints; j++)
04535                 checkvar[j] = -1;
04536         int k;
04537         for (int i = 0; i < numberOfStages; i++)
04538                 for (int j = 0; j < instanceData->timeDomain->stages->stage[i]->constraints->numberOfConstraints; j++)
04539                 {       k = instanceData->timeDomain->stages->stage[i]->constraints->con[j]->idx;
04540                         if (checkvar[k] != -1) 
04541                         {       delete [] checkvar;
04542                                 checkvar = NULL;
04543                                 return false;
04544                         }
04545                         checkvar[k] = instanceData->timeDomain->stages->stage[i]->constraints->con[j]->idx;
04546                 }
04547         delete [] checkvar;
04548         checkvar = NULL;                        
04549         return true;
04550 }// end setTimeDomainStageConstraintsUnordered()
04551 
04555 bool OSInstance::setTimeDomainStageObjectivesOrdered(int numberOfStages, int *numberOfObjectives, int *startIdx)
04556 {   if (instanceData->timeDomain == NULL)
04557                 instanceData->timeDomain = new TimeDomain();
04558         if (instanceData->timeDomain->interval != NULL)
04559                 return false;
04560         if (instanceData->timeDomain->stages == NULL)
04561                 instanceData->timeDomain->stages = new TimeDomainStages();      
04562         if (instanceData->timeDomain->stages != NULL)
04563         {       if ((instanceData->timeDomain->stages->numberOfStages != numberOfStages) &&
04564                         (instanceData->timeDomain->stages->numberOfStages != 0))
04565                          return false;
04566         }
04567         if (instanceData->timeDomain->stages->numberOfStages == 0)
04568         {       instanceData->timeDomain->stages->numberOfStages = numberOfStages;
04569                 if (instanceData->timeDomain->stages->stage == NULL)
04570                         instanceData->timeDomain->stages->stage = new TimeDomainStage*[numberOfStages];
04571                 for (int i = 0; i < numberOfStages; i++)
04572                         instanceData->timeDomain->stages->stage[i] = new TimeDomainStage();
04573         }
04574         for (int i = 0; i < numberOfStages; i++)
04575         //initial or empty vars, cons, objectives and set default to all objectives
04576         {       if (instanceData->timeDomain->stages->stage[i]->objectives != NULL)
04577                 {       delete instanceData->timeDomain->stages->stage[i]->objectives;
04578                         instanceData->timeDomain->stages->stage[i]->objectives = NULL;
04579                 }
04580                 instanceData->timeDomain->stages->stage[i]->objectives = new TimeDomainStageObjectives();
04581                 instanceData->timeDomain->stages->stage[i]->objectives->startIdx = startIdx[i];
04582                 instanceData->timeDomain->stages->stage[i]->objectives->numberOfObjectives = numberOfObjectives[i];
04583         }
04584         return true;
04585 }
04586 
04590 bool OSInstance::setTimeDomainStageObjectivesUnordered(int numberOfStages, int *numberOfObjectives, int **objIndex)
04591 {   if (instanceData->timeDomain == NULL)
04592                 instanceData->timeDomain = new TimeDomain();
04593         if (instanceData->timeDomain->interval != NULL)
04594                 return false;
04595         if (instanceData->timeDomain->stages == NULL)
04596                 instanceData->timeDomain->stages = new TimeDomainStages();      
04597         if (instanceData->timeDomain->stages != NULL)
04598         {       if ((instanceData->timeDomain->stages->numberOfStages != numberOfStages) &&
04599                         (instanceData->timeDomain->stages->numberOfStages != 0))
04600                          return false;
04601         }
04602         if (instanceData->timeDomain->stages->numberOfStages == 0)
04603         {       instanceData->timeDomain->stages->numberOfStages = numberOfStages;
04604                 if (instanceData->timeDomain->stages->stage == NULL)
04605                         instanceData->timeDomain->stages->stage = new TimeDomainStage*[numberOfStages];
04606                 for (int i = 0; i < numberOfStages; i++)
04607                         instanceData->timeDomain->stages->stage[i] = new TimeDomainStage();
04608         }
04609         for (int i = 0; i < numberOfStages; i++)
04610         //initial or empty vars, cons, objectives and set default to all objectives
04611         {       if (instanceData->timeDomain->stages->stage[i]->objectives != NULL)
04612                 {       delete instanceData->timeDomain->stages->stage[i]->objectives;
04613                         instanceData->timeDomain->stages->stage[i]->objectives = NULL;
04614                 }
04615                 instanceData->timeDomain->stages->stage[i]->objectives = new TimeDomainStageObjectives();
04616                 instanceData->timeDomain->stages->stage[i]->objectives->numberOfObjectives = numberOfObjectives[i];
04617                 instanceData->timeDomain->stages->stage[i]->objectives->obj = new TimeDomainStageObj*[numberOfObjectives[i]];
04618                 for (int j = 0; j < numberOfObjectives[i]; j++)
04619                 {       instanceData->timeDomain->stages->stage[i]->objectives->obj[j] = new TimeDomainStageObj();
04620                         instanceData->timeDomain->stages->stage[i]->objectives->obj[j]->idx = objIndex[i][j];
04621                 }
04622         }
04623         int *checkvar = new int[instanceData->objectives->numberOfObjectives];
04624         for (int j = 0; j < instanceData->objectives->numberOfObjectives; j++)
04625                 checkvar[j] = 0;
04626         int k;
04627         for (int i = 0; i < numberOfStages; i++)
04628                 for (int j = 0; j < instanceData->timeDomain->stages->stage[i]->objectives->numberOfObjectives; j++)
04629                 {       k = -instanceData->timeDomain->stages->stage[i]->objectives->obj[j]->idx-1;
04630                         checkvar[k] = instanceData->timeDomain->stages->stage[i]->objectives->obj[j]->idx;
04631                 }
04632         for (int i = 0; i < instanceData->objectives->numberOfObjectives; i++)
04633                 if (checkvar[i] == 0)
04634                 {       delete [] checkvar;
04635                         checkvar = NULL;
04636                         return false;
04637                 }
04638         delete [] checkvar;
04639         checkvar = NULL;                        
04640         return true;
04641 }
04642 
04646 bool OSInstance::setTimeDomainInterval(double start, double horizon)
04647 {   if (instanceData->timeDomain == NULL)
04648                 instanceData->timeDomain = new TimeDomain();
04649         if (instanceData->timeDomain->stages != NULL)
04650                 return false;
04651         if (instanceData->timeDomain->interval == NULL)
04652                 instanceData->timeDomain->interval = new TimeDomainInterval();  
04653         instanceData->timeDomain->interval->start = start;
04654         instanceData->timeDomain->interval->horizon = horizon;
04655         return true;
04656 } //end setTimeDomainInterval
04657 
04658 
04659 bool OSInstance::createOSADFun(std::vector<double> vdX){
04660         try{
04661                 if(m_bOSADFunIsCreated == true) return true;
04662                 //if( m_bNonLinearStructuresInitialized == false) initializeNonLinearStructures( );
04663                 if(m_binitForAlgDiff == false) initForAlgDiff();
04664                 
04665                 //if( m_bAllNonlinearVariablesIndex == false) getAllNonlinearVariablesIndexMap( );
04666                 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
04667                 unsigned int i;
04668                 size_t n = vdX.size();
04669 #ifdef COIN_HAS_CPPAD
04670                 // declare a CppAD vector and fill it in
04671                 CppAD::vector< AD<double> > vdaX( n );
04672                 for(i = 0; i < n; i++){
04673                         vdaX[ i] = vdX[ i];
04674                         //std::cout << "vdX =  " << vdX[ i] << std::endl;
04675                 }
04676                 // declare the independent variables and start recording
04677                 CppAD::Independent( vdaX);
04683                 CppAD::vector< AD<double> > m_vFG;        
04684                 int kount = 0;
04685                 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){  
04686                         m_vFG.push_back( (posMapExpTree->second)->m_treeRoot->constructADTape(&m_mapAllNonlinearVariablesIndex, &vdaX) );
04687                         //std::cout << "PUSHING BACK EXPRESSION NUMBER " << posMapExpTree->first << std::endl;
04688                         if( m_mapOSADFunRangeIndex.find( posMapExpTree->first) == m_mapOSADFunRangeIndex.end() ){
04689                                 // count which nonlinear obj/constraint this is
04690                                 m_mapOSADFunRangeIndex[ posMapExpTree->first] = kount;
04691                                 kount++;
04692                         }
04693                 }       
04694                 //create the function and stop recording
04695 //              std::cout << "create the function and stop recording"  << std::endl;
04696                 Fad = new CppAD::ADFun<double>(vdaX, m_vFG);
04697 //              std::cout << "range space dimension =  " << m_vFG.size() << std::endl;
04698                 // no forward sweeps done yet
04699                 m_iHighestTaylorCoeffOrder = -1;
04700                 m_bOSADFunIsCreated = true;
04701 #else
04702                 throw ErrorClass( "Error: An Algorithmic Differentiation Package Not Available");
04703 #endif
04704                         
04705                 return true;
04706         }
04707         catch(const ErrorClass& eclass){
04708                 throw ErrorClass( eclass.errormsg);
04709         } 
04710 }//end createOSADFun
04711 
04712 
04713 std::vector<double> OSInstance::forwardAD(int p, std::vector<double> vdX){
04714         try{
04715                 // make sure a OSADFun has been created
04716                 if(m_bOSADFunIsCreated == false) createOSADFun( vdX);
04717                 if(p > (m_iHighestTaylorCoeffOrder + 1) ) throw 
04718                         ErrorClass( "trying to calculate a p order forward when p-1 Taylor coefficient not available");
04719                 // adjust the order of the Taylor coefficient
04720                 m_iHighestTaylorCoeffOrder = p; 
04721                 m_iHighestOrderEvaluated = p;
04722                 //for(int i  = 0; i < vdX.size(); i++){
04723                         //std::cout << "ForwardAD Primal Variables " << i   << " " << vdX[ i] << std::endl;
04724                 //}
04725 #ifdef COIN_HAS_CPPAD
04726                 return (*Fad).Forward(p, vdX);
04727 #else
04728                 throw ErrorClass( "Error: An Algorithmic Differentiation Package Not Available");
04729 #endif
04730                 
04731         }
04732         catch(const ErrorClass& eclass){
04733                 throw ErrorClass( eclass.errormsg);
04734         }  
04735 }//end forwardAD
04736 
04737 
04738 std::vector<double> OSInstance::reverseAD(int p, std::vector<double> vdlambda){
04739         try{
04740 #ifndef COIN_HAS_CPPAD
04741                 throw ErrorClass( "Error: An Algorithmic Differentiation Package Not Available");
04742 #endif
04743                 if(p == 0) throw 
04744                         ErrorClass( "reverseAD must have p >= 1");
04745                 if(p > (m_iHighestTaylorCoeffOrder + 1) ) throw 
04746                         ErrorClass( "trying to calculate a p order reverse when p-1 Taylor coefficient not available");
04747                 //for(int i  = 0; i < vdlambda.size(); i++){
04748                 //      std::cout << "ReverseAD Multiplier " << i   << " " << vdlambda[ i] << std::endl;
04749                 //}
04750                 m_iHighestOrderEvaluated = p;
04751 #ifdef COIN_HAS_CPPAD
04752                 return (*Fad).Reverse(p, vdlambda);
04753 #else
04754                 throw ErrorClass( "Error: An Algorithmic Differentiation Package Not Available");
04755 #endif
04756                 
04757         }
04758         catch(const ErrorClass& eclass){
04759                 throw ErrorClass( eclass.errormsg);
04760         }  
04761 }//end forwardAD
04762 
04763 
04764 int  OSInstance::getADSparsityHessian(){
04765         
04766         unsigned int i;
04767         int numNonz;
04768         numNonz = 0;
04769         try{    
04770                 std::vector<bool> r(m_iNumberOfNonlinearVariables * m_iNumberOfNonlinearVariables);
04771                 unsigned int j;
04772                 for(i = 0; i < m_iNumberOfNonlinearVariables; i++) { 
04773                         for(j = 0; j < m_iNumberOfNonlinearVariables; j++)
04774                                 r[ i * m_iNumberOfNonlinearVariables + j ] = false;
04775                                 r[ i * m_iNumberOfNonlinearVariables + i] = true;
04776                 }       
04777                 // compute sparsity pattern for J(x) = F^{(1)} (x)
04778                 //should only be here if we have CppAD
04779 #ifdef COIN_HAS_CPPAD
04780                 (*Fad).ForSparseJac(m_iNumberOfNonlinearVariables, r);
04781 #else
04782                 throw ErrorClass( "Error: An Algorithmic Differentiation Package Not Available");
04783 #endif
04784                 //
04785                 //now the second derivative
04786                 unsigned int m = m_mapExpressionTreesMod.size();
04787                 std::vector<bool> e( m);
04788                 //Vector s(m);
04789                 for(i = 0; i < m; i++) e[i] = true;
04790 //              std::cout << "Computing Sparse Hessian" << std::endl;
04791                 //m_vbLagHessNonz holds the sparsity pattern Lagrangian of the Hessian
04792 #ifdef COIN_HAS_CPPAD
04793                 m_vbLagHessNonz = (*Fad).RevSparseHes(m_iNumberOfNonlinearVariables, e);
04794 #else
04795                 throw ErrorClass( "Error: An Algorithmic Differentiation Package Not Available");
04796 #endif
04797                 
04798                 for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
04799                         //std::cout << "Row " << i << "  of Hessian " << std::endl;
04800                         for(j = i; j < m_iNumberOfNonlinearVariables; j++){
04801                                 if(m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j]  == true) numNonz++;
04802                                 //std::cout << m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j] <<  "  " ;
04803                         }
04804                         //std::cout << std::endl;
04805                 }       
04806                 return numNonz;
04807         }
04808         catch(const ErrorClass& eclass){
04809                 throw ErrorClass( eclass.errormsg);
04810         }
04811 }//end getADSparsityHessian()
04812 
04813 

Generated on Mon Aug 3 03:02:22 2009 by  doxygen 1.4.7