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