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

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

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