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 timeDomain = NULL;
00756 }
00757
00758 string OSInstance::getInstanceName(){
00759 if( m_sInstanceName.length() <= 0){
00760 m_sInstanceName = instanceHeader->name;
00761 }
00762 return m_sInstanceName;
00763 }
00764
00765
00766 string OSInstance::getInstanceSource(){
00767 if( m_sInstanceSource.length() <= 0){
00768 m_sInstanceSource = instanceHeader->source;
00769 }
00770 return m_sInstanceSource;
00771 }
00772
00773 string OSInstance::getInstanceDescription(){
00774 if(m_sInstanceDescription.length() <= 0){
00775 m_sInstanceDescription = instanceHeader->description;
00776 }
00777 return m_sInstanceDescription;
00778 }
00779
00780 int OSInstance::getVariableNumber(){
00781 if(m_iVariableNumber == -1){
00782 m_iVariableNumber = instanceData->variables->numberOfVariables;
00783 }
00784 return m_iVariableNumber;
00785 }
00786
00787 int OSInstance::getNumberOfNonlinearExpressions(){
00788 if(m_iNonlinearExpressionNumber == -1){
00789 m_iNonlinearExpressionNumber = instanceData->nonlinearExpressions->numberOfNonlinearExpressions;
00790 }
00791 return m_iNonlinearExpressionNumber;
00792 }
00793
00794
00795
00796 bool OSInstance::processVariables() {
00797 if(m_bProcessVariables) return true;
00798 m_bProcessVariables = true;
00799 string vartype ="CBIS";
00800 int i = 0;
00801 int n = instanceData->variables->numberOfVariables;
00802 try{
00803 if(instanceData->variables->var == NULL) throw ErrorClass("no variables defined");
00804 if(instanceData->variables->var[0]->name.length() > 0 || instanceData->variables->var[n-1]->name.length() > 0){
00805 m_msVariableNames = new string[n];
00806 for(i = 0; i < n; i++) m_msVariableNames[i] = instanceData->variables->var[i]->name;
00807 }
00808 m_mdVariableInitialValues = new double[n];
00809 for(i = 0; i < n; i++){
00810 if(CommonUtil::ISOSNAN(instanceData->variables->var[ 0]->init) == true ){
00811 m_mdVariableInitialValues[i] = 1.7171;
00812 }
00813 else{
00814 m_mdVariableInitialValues[i] = instanceData->variables->var[i]->init;
00815 }
00816 }
00817
00818 if((instanceData->variables->var[0]->initString.length() > 0)){
00819 m_msVariableInitialStringValues = new string[n];
00820 for(i = 0; i < n; i++) m_msVariableInitialStringValues[i] = instanceData->variables->var[i]->initString;
00821 }
00822 m_mcVariableTypes = new char[n];
00823 m_mdVariableLowerBounds = new double[n];
00824 m_mdVariableUpperBounds = new double[n];
00825 for(i = 0; i < n; i++){
00826 if(vartype.find(instanceData->variables->var[i]->type) == string::npos) throw ErrorClass("wrong variable type");
00827 m_mcVariableTypes[i] = instanceData->variables->var[i]->type;
00828 if(m_mcVariableTypes[i] == 'B') m_iNumberOfBinaryVariables++;
00829 if(m_mcVariableTypes[i] == 'I') m_iNumberOfIntegerVariables++;
00830 m_mdVariableLowerBounds[i] = instanceData->variables->var[i]->lb;
00831 m_mdVariableUpperBounds[i] = instanceData->variables->var[i]->ub;
00832 }
00833 return true;
00834 }
00835 catch(const ErrorClass& eclass){
00836 throw ErrorClass( eclass.errormsg);
00837 }
00838 }
00839
00840 string* OSInstance::getVariableNames() {
00841 processVariables();
00842 return m_msVariableNames;
00843 }
00844
00845 double* OSInstance::getVariableInitialValues() {
00846 processVariables();
00847 return m_mdVariableInitialValues;
00848 }
00849
00850 string* OSInstance::getVariableInitialStringValues() {
00851 processVariables();
00852 return m_msVariableInitialStringValues;
00853 }
00854
00855 char* OSInstance::getVariableTypes() {
00856 processVariables();
00857 return m_mcVariableTypes;
00858 }
00859
00860 int OSInstance::getNumberOfIntegerVariables() {
00861 processVariables();
00862 return m_iNumberOfIntegerVariables;
00863 }
00864
00865 int OSInstance::getNumberOfBinaryVariables() {
00866 processVariables();
00867 return m_iNumberOfBinaryVariables;
00868 }
00869
00870 double* OSInstance::getVariableLowerBounds() {
00871 processVariables();
00872 return m_mdVariableLowerBounds;
00873 }
00874
00875 double* OSInstance::getVariableUpperBounds() {
00876 processVariables();
00877 return m_mdVariableUpperBounds;
00878 }
00879
00880 int OSInstance::getObjectiveNumber(){
00881 if(m_iObjectiveNumber == -1){
00882 m_iObjectiveNumber = instanceData->objectives->numberOfObjectives;
00883 }
00884 return m_iObjectiveNumber;
00885 }
00886
00887
00888 bool OSInstance::processObjectives() {
00889 if(m_bProcessObjectives) return true;
00890 m_bProcessObjectives = true;
00891 int i = 0;
00892 int j = 0;
00893 if(instanceData == NULL || instanceData->objectives == NULL || instanceData->objectives->obj == NULL || instanceData->objectives->numberOfObjectives == 0) return true;
00894 int n = instanceData->objectives->numberOfObjectives;
00895 try{
00896 if(instanceData->objectives->obj[0]->name.length() > 0 || instanceData->objectives->obj[n-1]->name.length() > 0){
00897 m_msObjectiveNames = new string[n];
00898 for(i = 0; i < n; i++) m_msObjectiveNames[i] = instanceData->objectives->obj[i]->name;
00899 }
00900 m_msMaxOrMins = new string[n];
00901 m_miNumberOfObjCoef = new int[n];
00902 m_mdObjectiveConstants = new double[n];
00903 m_mdObjectiveWeights = new double[n];
00904 m_mObjectiveCoefficients = new SparseVector*[n];
00905 for(i = 0; i < n; i++){
00906 m_mObjectiveCoefficients[i] = new SparseVector(instanceData->objectives->obj[ j]->numberOfObjCoef);
00907
00908 }
00909
00910
00911
00912
00913
00914 for(i = 0; i < n; i++){
00915 if((instanceData->objectives->obj[i]->maxOrMin.compare("max") != 0) && (instanceData->objectives->obj[i]->maxOrMin.compare("min") != 0 )) throw ErrorClass("wrong objective maxOrMin");
00916 m_msMaxOrMins[i] = instanceData->objectives->obj[i]->maxOrMin;
00917 m_miNumberOfObjCoef[i] = instanceData->objectives->obj[i]->numberOfObjCoef;
00918 m_mdObjectiveConstants[i] = instanceData->objectives->obj[i]->constant;
00919 m_mdObjectiveWeights[i] = instanceData->objectives->obj[i]->weight;
00920 if(instanceData->objectives->obj[i]->coef == NULL && m_miNumberOfObjCoef[i] != 0){
00921 throw ErrorClass("objective coefficient number inconsistent with objective coefficient array");
00922 }
00923 if(instanceData->objectives->obj[i]->coef != NULL){
00924 for(j = 0; j < m_mObjectiveCoefficients[i]->number; j++){
00925 m_mObjectiveCoefficients[i]->indexes[j] = instanceData->objectives->obj[i]->coef[j]->idx;
00926 m_mObjectiveCoefficients[i]->values[j] = instanceData->objectives->obj[i]->coef[j]->value;
00927 }
00928 }
00929 }
00930 return true;
00931 }
00932 catch(const ErrorClass& eclass){
00933 throw ErrorClass( eclass.errormsg);
00934 }
00935 }
00936
00937 string* OSInstance::getObjectiveNames() {
00938 processObjectives();
00939 return m_msObjectiveNames;
00940 }
00941
00942 string* OSInstance::getObjectiveMaxOrMins() {
00943 processObjectives();
00944 return m_msMaxOrMins;
00945 }
00946
00947 int* OSInstance::getObjectiveCoefficientNumbers(){
00948 processObjectives();
00949 return m_miNumberOfObjCoef;
00950 }
00951
00952 double* OSInstance::getObjectiveConstants() {
00953 processObjectives();
00954 return m_mdObjectiveConstants;
00955 }
00956
00957 double* OSInstance::getObjectiveWeights() {
00958 processObjectives();
00959 return m_mdObjectiveWeights;
00960 }
00961
00962 SparseVector** OSInstance::getObjectiveCoefficients() {
00963 processObjectives();
00964 return m_mObjectiveCoefficients;
00965 }
00966
00967
00968 double** OSInstance::getDenseObjectiveCoefficients() {
00969 if(m_bGetDenseObjectives) return m_mmdDenseObjectiveCoefficients;
00970 m_bGetDenseObjectives = true;
00971 if(instanceData->objectives->obj == NULL || instanceData->objectives->numberOfObjectives == 0) return m_mmdDenseObjectiveCoefficients;
00972 int m = instanceData->objectives->numberOfObjectives;
00973 int n = instanceData->variables->numberOfVariables;
00974 m_mmdDenseObjectiveCoefficients = new double*[m];
00975 int i, j, numobjcoef;
00976 SparseVector *sparsevec;
00977 for(i = 0; i < m; i++){
00978 sparsevec = this->getObjectiveCoefficients()[i];
00979 m_mmdDenseObjectiveCoefficients[ i] = new double[n];
00980 for(j = 0; j < n; j++){
00981 m_mmdDenseObjectiveCoefficients[ i][j] = 0.0;
00982 }
00983 sparsevec = this->getObjectiveCoefficients()[i];
00984 numobjcoef = sparsevec->number;
00985 for(j = 0; j < numobjcoef; j++){
00986 m_mmdDenseObjectiveCoefficients[i][ sparsevec->indexes[ j]]
00987 = sparsevec->values[ j];
00988 }
00989 }
00990 return m_mmdDenseObjectiveCoefficients;
00991 }
00992
00993 int OSInstance::getConstraintNumber(){
00994 if(m_iConstraintNumber == -1){
00995 m_iConstraintNumber = instanceData->constraints->numberOfConstraints;
00996 }
00997 return m_iConstraintNumber;
00998 }
00999
01000 bool OSInstance::processConstraints() {
01001 if(m_bProcessConstraints) return true;
01002 m_bProcessConstraints = true;
01003 int i = 0;
01004 ostringstream outStr;
01005 if(instanceData == NULL || instanceData->constraints == NULL || instanceData->constraints->con == NULL || instanceData->constraints->numberOfConstraints == 0) return true;
01006 int n = instanceData->constraints->numberOfConstraints;
01007 try{
01008 if(instanceData->constraints->con[0]->name.length() > 0 || instanceData->constraints->con[n-1]->name.length() > 0){
01009 m_msConstraintNames = new string[n];
01010 for(i = 0; i < n; i++) m_msConstraintNames[i] = instanceData->constraints->con[i]->name;
01011 }
01012 m_mdConstraintLowerBounds = new double[n];
01013 m_mdConstraintUpperBounds = new double[n];
01014 m_mdConstraintConstants = new double[n];
01015 m_mcConstraintTypes = new char[n];
01016 for(i = 0; i < n; i++){
01017 m_mdConstraintLowerBounds[i] = instanceData->constraints->con[i]->lb;
01018 m_mdConstraintUpperBounds[i] = instanceData->constraints->con[i]->ub;
01019 m_mdConstraintConstants[i] = instanceData->constraints->con[i]->constant;
01020 if(m_mdConstraintLowerBounds[i] == OSDBL_MAX || m_mdConstraintUpperBounds[i] == -OSDBL_MAX) {
01021 throw ErrorClass( outStr.str() );
01022 }
01023 else if(m_mdConstraintLowerBounds[i] > m_mdConstraintUpperBounds[i]) {
01024 outStr << "Constraint " ;
01025 outStr << i;
01026 outStr << " is infeasible";
01027 throw ErrorClass( outStr.str());
01028 }
01029 else if(m_mdConstraintLowerBounds[i] == -OSDBL_MAX && m_mdConstraintUpperBounds[i] == OSDBL_MAX)
01030 m_mcConstraintTypes[i] = 'U';
01031 else if(m_mdConstraintLowerBounds[i] == m_mdConstraintUpperBounds[i])
01032 m_mcConstraintTypes[i] = 'E';
01033 else if(m_mdConstraintLowerBounds[i] == -OSDBL_MAX)
01034 m_mcConstraintTypes[i] = 'L';
01035 else if(m_mdConstraintUpperBounds[i] == OSDBL_MAX)
01036 m_mcConstraintTypes[i] = 'G';
01037 else m_mcConstraintTypes[i] = 'R';
01038 }
01039 return true;
01040 }
01041 catch(const ErrorClass& eclass){
01042 throw ErrorClass( eclass.errormsg);
01043 }
01044 }
01045
01046
01047 string* OSInstance::getConstraintNames() {
01048 processConstraints();
01049 return m_msConstraintNames;
01050 }
01051
01052
01053 double* OSInstance::getConstraintLowerBounds() {
01054 processConstraints();
01055 return m_mdConstraintLowerBounds;
01056 }
01057
01058 char* OSInstance::getConstraintTypes() {
01059 processConstraints();
01060 return m_mcConstraintTypes;
01061 }
01062
01063 double* OSInstance::getConstraintUpperBounds() {
01064 processConstraints();
01065 return m_mdConstraintUpperBounds;
01066 }
01067
01068 int OSInstance::getLinearConstraintCoefficientNumber(){
01069 if(m_iLinearConstraintCoefficientNumber == -1){
01070 m_iLinearConstraintCoefficientNumber = instanceData->linearConstraintCoefficients->numberOfValues;
01071 }
01072 return m_iLinearConstraintCoefficientNumber;
01073 }
01074
01075 bool OSInstance::processLinearConstraintCoefficients() {
01076 if(m_bProcessLinearConstraintCoefficients) return true;
01077 m_bProcessLinearConstraintCoefficients = true;
01078 try{
01079 int n = instanceData->linearConstraintCoefficients->numberOfValues;
01080
01081 if((instanceData->linearConstraintCoefficients->value == NULL ) || (n == 0) ) return true;
01082
01083 if((instanceData->linearConstraintCoefficients->colIdx != NULL && instanceData->linearConstraintCoefficients->colIdx->el != NULL)
01084 && (instanceData->linearConstraintCoefficients->rowIdx != NULL && instanceData->linearConstraintCoefficients->rowIdx->el != NULL))
01085 throw ErrorClass("ambiguous linear constraint coefficient major");
01086 else if(instanceData->linearConstraintCoefficients->value->el == NULL) return true;
01087 else{
01088 if(instanceData->linearConstraintCoefficients->rowIdx->el != NULL){
01089 m_bColumnMajor = true;
01090 m_linearConstraintCoefficientsInColumnMajor = new SparseMatrix();
01091 m_linearConstraintCoefficientsInColumnMajor->bDeleteArrays = false;
01092 m_linearConstraintCoefficientsInColumnMajor->isColumnMajor = true;
01093 m_linearConstraintCoefficientsInColumnMajor->valueSize = n;
01094 m_linearConstraintCoefficientsInColumnMajor->startSize = instanceData->variables->numberOfVariables + 1;
01095 }
01096 else{
01097 m_bColumnMajor = false;
01098 m_linearConstraintCoefficientsInRowMajor = new SparseMatrix();
01099 m_linearConstraintCoefficientsInRowMajor->bDeleteArrays = false;
01100 m_linearConstraintCoefficientsInRowMajor->isColumnMajor = false;
01101 m_linearConstraintCoefficientsInRowMajor->valueSize = n;
01102 m_linearConstraintCoefficientsInRowMajor->startSize = instanceData->constraints->numberOfConstraints + 1;
01103 }
01104 }
01105 if(m_bColumnMajor == true){
01106 m_linearConstraintCoefficientsInColumnMajor->values = instanceData->linearConstraintCoefficients->value->el;
01107 m_linearConstraintCoefficientsInColumnMajor->indexes = instanceData->linearConstraintCoefficients->rowIdx->el;
01108 m_linearConstraintCoefficientsInColumnMajor->starts = instanceData->linearConstraintCoefficients->start->el;
01109 }
01110 else{
01111 m_linearConstraintCoefficientsInRowMajor->values = instanceData->linearConstraintCoefficients->value->el;
01112 m_linearConstraintCoefficientsInRowMajor->indexes = instanceData->linearConstraintCoefficients->colIdx->el;
01113 m_linearConstraintCoefficientsInRowMajor->starts = instanceData->linearConstraintCoefficients->start->el;
01114 }
01115 return true;
01116 }
01117 catch(const ErrorClass& eclass){
01118 throw ErrorClass( eclass.errormsg);
01119 }
01120 }
01121
01122 bool OSInstance::getLinearConstraintCoefficientMajor() {
01123 processLinearConstraintCoefficients();
01124 return m_bColumnMajor;
01125 }
01126
01127 SparseMatrix* OSInstance::getLinearConstraintCoefficientsInColumnMajor() {
01128 processLinearConstraintCoefficients();
01129 if(m_linearConstraintCoefficientsInColumnMajor != NULL) return m_linearConstraintCoefficientsInColumnMajor;
01130 if(!m_bColumnMajor){
01131 if(m_linearConstraintCoefficientsInRowMajor == NULL) return NULL;
01132 m_linearConstraintCoefficientsInColumnMajor =
01133 MathUtil::convertLinearConstraintCoefficientMatrixToTheOtherMajor(false,
01134 m_linearConstraintCoefficientsInRowMajor->startSize,
01135 m_linearConstraintCoefficientsInRowMajor->valueSize,
01136 m_linearConstraintCoefficientsInRowMajor->starts,
01137 m_linearConstraintCoefficientsInRowMajor->indexes,
01138 m_linearConstraintCoefficientsInRowMajor->values,
01139 getVariableNumber());
01140 }
01141 return m_linearConstraintCoefficientsInColumnMajor;
01142 }
01143
01144 SparseMatrix* OSInstance::getLinearConstraintCoefficientsInRowMajor() {
01145 processLinearConstraintCoefficients();
01146 if(m_linearConstraintCoefficientsInRowMajor != NULL) return m_linearConstraintCoefficientsInRowMajor;
01147 if(m_bColumnMajor){
01148 if(m_linearConstraintCoefficientsInColumnMajor == NULL) return NULL;
01149 m_linearConstraintCoefficientsInRowMajor =
01150 MathUtil::convertLinearConstraintCoefficientMatrixToTheOtherMajor(true,
01151 m_linearConstraintCoefficientsInColumnMajor->startSize,
01152 m_linearConstraintCoefficientsInColumnMajor->valueSize,
01153 m_linearConstraintCoefficientsInColumnMajor->starts,
01154 m_linearConstraintCoefficientsInColumnMajor->indexes,
01155 m_linearConstraintCoefficientsInColumnMajor->values,
01156 getConstraintNumber());
01157 }
01158 return m_linearConstraintCoefficientsInRowMajor;
01159 }
01160
01161
01162 int OSInstance::getNumberOfQuadraticTerms(){
01163 if(m_iQuadraticTermNumber == -1){
01164
01165
01166 if(instanceData->quadraticCoefficients == NULL)instanceData->quadraticCoefficients = new QuadraticCoefficients();
01167 m_iQuadraticTermNumber = instanceData->quadraticCoefficients->numberOfQuadraticTerms;
01168 }
01169 return m_iQuadraticTermNumber;
01170 }
01171
01172 QuadraticTerms* OSInstance::getQuadraticTerms() {
01173 if(m_bProcessQuadraticTerms) return m_quadraticTerms;
01174 m_bProcessQuadraticTerms = true;
01175 if(instanceData->quadraticCoefficients->qTerm == 0) return 0;
01176 try{
01177 int i = 0;
01178 QuadraticCoefficients* quadraticCoefs = instanceData->quadraticCoefficients;
01179 int n = quadraticCoefs->numberOfQuadraticTerms;
01180 if(!quadraticCoefs->qTerm && n != 0)
01181 throw ErrorClass("quadratic term number inconsistent with quadratic term array");
01182 m_quadraticTerms = new QuadraticTerms();
01183 m_quadraticTerms->rowIndexes = new int[n];
01184 m_quadraticTerms->varOneIndexes = new int[n];
01185 m_quadraticTerms->varTwoIndexes = new int[n];
01186 m_quadraticTerms->coefficients = new double[n];
01187 for(i = 0; i < n; i++){
01188 m_quadraticTerms->rowIndexes[i] = quadraticCoefs->qTerm[i]->idx;
01189 m_quadraticTerms->varOneIndexes[i] = quadraticCoefs->qTerm[i]->idxOne;
01190 m_quadraticTerms->varTwoIndexes[i] = quadraticCoefs->qTerm[i]->idxTwo;
01191 m_quadraticTerms->coefficients[i] = quadraticCoefs->qTerm[i]->coef;
01192 }
01193 return m_quadraticTerms;
01194 }
01195 catch(const ErrorClass& eclass){
01196 throw ErrorClass( eclass.errormsg);
01197 }
01198 }
01199
01200
01201 int* OSInstance::getQuadraticRowIndexes() {
01202 if(m_bQuadraticRowIndexesProcessed == true) return m_miQuadRowIndexes;
01203 m_bQuadraticRowIndexesProcessed = true;
01204 int n = getNumberOfQuadraticTerms();
01205 if(n <= 0) return NULL;
01206 QuadraticTerms *qTerms = NULL;
01207 qTerms = getQuadraticTerms();
01208 std::map<int, int> foundIdx;
01209 std::map<int, int>::iterator pos;
01210 int i;
01211 try{
01212 for(i = 0; i < n; i++){
01213
01214 foundIdx[ qTerms->rowIndexes[ i] ];
01215 }
01216
01217 m_iNumberOfQuadraticRowIndexes = foundIdx.size();
01218 m_miQuadRowIndexes = new int[ m_iNumberOfQuadraticRowIndexes ] ;
01219 i = 0;
01220 for(pos = foundIdx.begin(); pos != foundIdx.end(); ++pos){
01221 m_miQuadRowIndexes[ i++] = pos->first;
01222 }
01223 foundIdx.clear();
01224 return m_miQuadRowIndexes;
01225 }
01226 catch(const ErrorClass& eclass){
01227 throw ErrorClass( eclass.errormsg);
01228 }
01229 }
01230
01231
01232 int OSInstance::getNumberOfQuadraticRowIndexes() {
01233 if(m_bQuadraticRowIndexesProcessed == false) getQuadraticRowIndexes();
01234 return m_iNumberOfQuadraticRowIndexes;
01235 }
01236
01237 int* OSInstance::getNonlinearExpressionTreeIndexes(){
01238 if(m_bNonlinearExpressionTreeIndexesProcessed == true) return m_miNonlinearExpressionTreeIndexes;
01239 m_bNonlinearExpressionTreeIndexesProcessed = true;
01240 std::map<int, OSExpressionTree*> expTrees;
01241 expTrees = getAllNonlinearExpressionTrees();
01242 std::map<int, OSExpressionTree*>::iterator pos;
01243 try{
01244
01245 m_iNumberOfNonlinearExpressionTreeIndexes = expTrees.size();
01246 m_miNonlinearExpressionTreeIndexes = new int[ m_iNumberOfNonlinearExpressionTreeIndexes ] ;
01247 int i = 0;
01248 for(pos = expTrees.begin(); pos != expTrees.end(); ++pos){
01249 m_miNonlinearExpressionTreeIndexes[ i++] = pos->first;
01250 }
01251 expTrees.clear();
01252 return m_miNonlinearExpressionTreeIndexes;
01253 }
01254 catch(const ErrorClass& eclass){
01255 throw ErrorClass( eclass.errormsg);
01256 }
01257 }
01258
01259 int OSInstance::getNumberOfNonlinearExpressionTreeIndexes() {
01260 if(m_bNonlinearExpressionTreeIndexesProcessed == false) getNonlinearExpressionTreeIndexes();
01261 return m_iNumberOfNonlinearExpressionTreeIndexes;
01262 }
01263
01264
01265
01266 int* OSInstance::getNonlinearExpressionTreeModIndexes(){
01267 if(m_bNonlinearExpressionTreeModIndexesProcessed == true) return m_miNonlinearExpressionTreeModIndexes;
01268 m_bNonlinearExpressionTreeModIndexesProcessed = true;
01269 std::map<int, OSExpressionTree*> expTrees;
01270 expTrees = getAllNonlinearExpressionTreesMod();
01271 std::map<int, OSExpressionTree*>::iterator pos;
01272 try{
01273
01274 m_iNumberOfNonlinearExpressionTreeModIndexes = expTrees.size();
01275 m_miNonlinearExpressionTreeModIndexes = new int[ m_iNumberOfNonlinearExpressionTreeModIndexes ] ;
01276 int i = 0;
01277 for(pos = expTrees.begin(); pos != expTrees.end(); ++pos){
01278 m_miNonlinearExpressionTreeModIndexes[ i++] = pos->first;
01279 }
01280 expTrees.clear();
01281 return m_miNonlinearExpressionTreeModIndexes;
01282 }
01283 catch(const ErrorClass& eclass){
01284 throw ErrorClass( eclass.errormsg);
01285 }
01286 }
01287
01288 int OSInstance::getNumberOfNonlinearExpressionTreeModIndexes() {
01289 if(m_bNonlinearExpressionTreeModIndexesProcessed == false) getNonlinearExpressionTreeModIndexes();
01290 return m_iNumberOfNonlinearExpressionTreeModIndexes;
01291 }
01292
01293
01294 int OSInstance::getNumberOfNonlinearConstraints(){
01295 if( m_bProcessExpressionTrees == false ){
01296 getAllNonlinearExpressionTrees();
01297 return m_iConstraintNumberNonlinear;
01298 }
01299 else return m_iConstraintNumberNonlinear;
01300 }
01301
01302 int OSInstance::getNumberOfNonlinearObjectives(){
01303 if( m_bProcessExpressionTrees == false ){
01304 getAllNonlinearExpressionTrees();
01305 return m_iObjectiveNumberNonlinear;
01306 }
01307 else return m_iObjectiveNumberNonlinear;
01308 }
01309
01310
01311 OSExpressionTree* OSInstance::getNonlinearExpressionTree(int rowIdx){
01312 if( m_bProcessExpressionTrees == false ){
01313 getAllNonlinearExpressionTrees();
01314 return m_mapExpressionTrees[ rowIdx];
01315 }
01316 else{
01317
01318
01319 if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()) return m_mapExpressionTrees[ rowIdx];
01320 else return NULL ;
01321
01323
01325
01326
01327
01328
01329 }
01330 }
01331
01332
01333 OSExpressionTree* OSInstance::getNonlinearExpressionTreeMod(int rowIdx){
01334 if( m_bProcessExpressionTreesMod == false ){
01335 getAllNonlinearExpressionTreesMod();
01336 return m_mapExpressionTreesMod[ rowIdx];
01337 }
01338 else{
01339
01340
01341 if( m_mapExpressionTreesMod.find( rowIdx) != m_mapExpressionTreesMod.end()) return m_mapExpressionTreesMod[ rowIdx];
01342 else return NULL ;
01343
01345
01347
01348
01349
01350
01351 }
01352 }
01353
01354
01355 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeInPostfix( int rowIdx){
01356 if( m_bProcessExpressionTrees == false ) getAllNonlinearExpressionTrees();
01357 std::vector<OSnLNode*> postfixVec;
01358 try{
01359 if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()){
01360 OSExpressionTree* expTree = getNonlinearExpressionTree( rowIdx);
01361 postfixVec = expTree->m_treeRoot->getPostfixFromExpressionTree();
01362
01363 }
01364 else{
01365 throw ErrorClass("Error in getNonlinearExpressionTreeInPostfix, rowIdx not valid");
01366 }
01367 return postfixVec;
01368 }
01369 catch(const ErrorClass& eclass){
01370 throw ErrorClass( eclass.errormsg);
01371 }
01372 }
01373
01374
01375 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeModInPostfix( int rowIdx){
01376 if( m_bProcessExpressionTreesMod == false ) getAllNonlinearExpressionTreesMod();
01377 std::vector<OSnLNode*> postfixVec;
01378 try{
01379 if( m_mapExpressionTreesMod.find( rowIdx) != m_mapExpressionTreesMod.end()){
01380 OSExpressionTree* expTree = getNonlinearExpressionTreeMod( rowIdx);
01381 postfixVec = expTree->m_treeRoot->getPostfixFromExpressionTree();
01382
01383 }
01384 else{
01385 throw ErrorClass("Error in getNonlinearExpressionTreeModInPostfix, rowIdx not valid");
01386 }
01387 return postfixVec;
01388 }
01389 catch(const ErrorClass& eclass){
01390 throw ErrorClass( eclass.errormsg);
01391 }
01392 }
01393
01394
01395 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeInPrefix( int rowIdx){
01396 if( m_bProcessExpressionTrees == false ) getAllNonlinearExpressionTrees();
01397 std::vector<OSnLNode*> prefixVec;
01398 try{
01399 if( m_mapExpressionTrees.find( rowIdx) != m_mapExpressionTrees.end()){
01400 OSExpressionTree* expTree = getNonlinearExpressionTree( rowIdx);
01401 prefixVec = expTree->m_treeRoot->getPrefixFromExpressionTree();
01402
01403 }
01404 else{
01405 throw ErrorClass("Error in getNonlinearExpressionTreeInPrefix, rowIdx not valid");
01406 }
01407 return prefixVec;
01408 }
01409 catch(const ErrorClass& eclass){
01410 throw ErrorClass( eclass.errormsg);
01411 }
01412 }
01413
01414 std::vector<OSnLNode*> OSInstance::getNonlinearExpressionTreeModInPrefix( int rowIdx){
01415 if( m_bProcessExpressionTreesMod == false ) getAllNonlinearExpressionTreesMod();
01416 std::vector<OSnLNode*> prefixVec;
01417 try{
01418 if( m_mapExpressionTreesMod.find( rowIdx) != m_mapExpressionTreesMod.end()){
01419 OSExpressionTree* expTree = getNonlinearExpressionTreeMod( rowIdx);
01420 prefixVec = expTree->m_treeRoot->getPrefixFromExpressionTree();
01421
01422 }
01423 else{
01424 throw ErrorClass("Error in getNonlinearExpressionTreeInPrefix, rowIdx not valid");
01425 }
01426 return prefixVec;
01427 }
01428 catch(const ErrorClass& eclass){
01429 throw ErrorClass( eclass.errormsg);
01430 }
01431 }
01432
01433 std::map<int, OSExpressionTree*> OSInstance::getAllNonlinearExpressionTrees(){
01434 if(m_bProcessExpressionTrees == true) return m_mapExpressionTrees;
01435 std::map<int, int> foundIdx;
01436 std::map<int, int>::iterator pos;
01437 OSnLNodePlus *nlNodePlus;
01438 OSExpressionTree *expTree;
01439 m_iObjectiveNumberNonlinear = 0;
01440 m_iConstraintNumberNonlinear = 0;
01441 int i;
01442
01443 if( instanceData->nonlinearExpressions->numberOfNonlinearExpressions > 0 && instanceData->nonlinearExpressions->nl != NULL){
01444 for( i = 0; i < instanceData->nonlinearExpressions->numberOfNonlinearExpressions; i++){
01445 instanceData->nonlinearExpressions->nl[i]->m_bDeleteExpressionTree = false;
01446 }
01447 }
01448 int index;
01449
01450 for(i = 0; i < instanceData->nonlinearExpressions->numberOfNonlinearExpressions; i++){
01451 index = instanceData->nonlinearExpressions->nl[ i]->idx;
01452 if(foundIdx.find( index) != foundIdx.end() ){
01453
01454
01455
01456
01457
01458
01459
01460 nlNodePlus = new OSnLNodePlus();
01461
01462 expTree = instanceData->nonlinearExpressions->nl[ i]->osExpressionTree;
01463
01464 nlNodePlus->m_mChildren[ 0] = m_mapExpressionTrees[ index]->m_treeRoot;
01465 nlNodePlus->m_mChildren[ 1] = instanceData->nonlinearExpressions->nl[ i]->osExpressionTree->m_treeRoot;
01466
01467 instanceData->nonlinearExpressions->nl[ foundIdx[ index] ]->m_bDeleteExpressionTree = true;
01468 instanceData->nonlinearExpressions->nl[ foundIdx[ index] ]->osExpressionTree->bDestroyNlNodes = false;
01469
01470 m_mapExpressionTrees[ index] = expTree;
01471 m_mapExpressionTrees[ index]->m_treeRoot = nlNodePlus;
01472 foundIdx[ index] = i;
01473 }
01474 else{
01475
01476 m_mapExpressionTrees[ index] = instanceData->nonlinearExpressions->nl[ i]->osExpressionTree;
01477 m_mapExpressionTrees[ index]->m_treeRoot = instanceData->nonlinearExpressions->nl[ i]->osExpressionTree->m_treeRoot;
01478 foundIdx[ index] = i;
01479 }
01480
01481 }
01482
01483 for(pos = foundIdx.begin(); pos != foundIdx.end(); ++pos){
01484 if(pos->first == -1) {
01485 m_iObjectiveNumberNonlinear++;
01486 }
01487 else m_iConstraintNumberNonlinear++;
01488 }
01489 foundIdx.clear();
01490 m_bProcessExpressionTrees = true;
01491 return m_mapExpressionTrees;
01492 }
01493
01494 std::map<int, OSExpressionTree*> OSInstance::getAllNonlinearExpressionTreesMod(){
01495 if( m_bProcessExpressionTreesMod == true ) return m_mapExpressionTreesMod;
01496 m_bProcessExpressionTreesMod = true;
01497
01498 if( m_bNonLinearStructuresInitialized == false) initializeNonLinearStructures( );
01499 return m_mapExpressionTreesMod;
01500 }
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512 bool OSInstance::setInstanceSource(string source){
01513 instanceHeader->source = source;
01514 return true;
01515 }
01516
01517 bool OSInstance::setInstanceDescription(string description){
01518 instanceHeader->description = description;
01519 return true;
01520 }
01521
01522
01523 bool OSInstance::setInstanceName(string name){
01524 instanceHeader->name = name;
01525 return true;
01526 }
01527
01528
01529 bool OSInstance::setVariableNumber(int number){
01530
01531 if(number <= 0) return false;
01532
01533
01534
01535
01536 if(instanceData->variables == NULL) instanceData->variables = new Variables();
01537 instanceData->variables->numberOfVariables = number;
01538 if(instanceData->variables->var == NULL){
01539 instanceData->variables->var = new Variable*[number];
01540 }
01541 return true;
01542 }
01543
01544
01545 bool OSInstance::addVariable(int index, string name, double lowerBound, double upperBound, char type, double init, string initString){
01546 instanceData->variables->var[index] = new Variable();
01547 if(index < 0 || instanceData->variables->numberOfVariables <= 0 || index >= instanceData->variables->numberOfVariables) return false;
01548 instanceData->variables->var[index]->name = name;
01549
01550 instanceData->variables->var[index]->lb = lowerBound;
01551 if(upperBound != OSDBL_MAX && upperBound != OSDBL_MAX)instanceData->variables->var[index]->ub = upperBound;
01552 instanceData->variables->var[index]->type = type;
01553 if(init != OSNAN) instanceData->variables->var[index]->init = init;
01554 instanceData->variables->var[index]->initString = initString;
01555 return true;
01556 }
01557
01558
01559 bool OSInstance::setVariables(int number, string *names, double *lowerBounds,
01560 double *upperBounds, char *types, double *inits, string *initsString){
01561 if(number <= 0) return false;
01562 try{
01563 if(instanceData->variables == NULL){
01564 throw ErrorClass("There is no variables object");
01565 }
01566 if(instanceData->variables->numberOfVariables != number){
01567 throw ErrorClass("input number of variables not equal to number in class");
01568 }
01569
01570 int i;
01571 for(i = 0; i < number; i++){
01572 instanceData->variables->var[ i] = new Variable();
01573 }
01574 if(names != NULL){
01575 for(i = 0; i < number; i++) instanceData->variables->var[i]->name = names[i];
01576 }
01577 if(lowerBounds != NULL){
01578 for(i = 0; i < number; i++){
01579 if(lowerBounds[i] != -OSDBL_MAX && lowerBounds[i] != -OSDBL_MAX)instanceData->variables->var[i]->lb = lowerBounds[i];
01580 }
01581 }
01582 if(upperBounds != NULL){
01583 for(i = 0; i < number; i++){
01584 if(upperBounds[i] != OSDBL_MAX && upperBounds[i] != OSDBL_MAX)instanceData->variables->var[i]->ub = upperBounds[i];
01585 }
01586 }
01587 if(types != NULL){
01588 for(i = 0; i < number; i++){
01589 instanceData->variables->var[i]->type = types[i];
01590 if(types[i] != 'C' && types[i] != 'B' && types[i] != 'I' && types[i] != 'S') types[i] = 'C';
01591 }
01592 }
01593 if(inits != NULL){
01594 for(i = 0; i < number; i++) instanceData->variables->var[i]->init = inits[i];
01595 }
01596 if(initsString != NULL){
01597 for(i = 0; i < number; i++) instanceData->variables->var[i]->initString = initsString[i];
01598 }
01599 return true;
01600 }
01601 catch(const ErrorClass& eclass){
01602 throw ErrorClass( eclass.errormsg);
01603 }
01604 }
01605
01606
01607
01608 bool OSInstance::setObjectiveNumber(int number){
01609 if(number < 0) return false;
01610 if(instanceData->objectives == NULL) instanceData->objectives = new Objectives();
01611 if(number == 0){
01612 instanceData->objectives->numberOfObjectives = 0;
01613 instanceData->objectives->obj = 0;
01614 return true;
01615 }
01616 instanceData->objectives->numberOfObjectives = number;
01617 instanceData->objectives->obj = new Objective*[number];
01618 return true;
01619 }
01620
01621 bool OSInstance::addObjective(int index, string name, string maxOrMin, double constant, double weight, SparseVector *objectiveCoefficients){
01622 if(index >= 0 || instanceData->objectives->numberOfObjectives <= 0 || abs(index) > instanceData->objectives->numberOfObjectives) return false;
01623 int arrayIndex = abs(index) -1;
01624 if(instanceData->objectives->obj == NULL) return false;
01625 instanceData->objectives->obj[arrayIndex] = new Objective();
01626 instanceData->objectives->obj[arrayIndex]->name = name;
01627 if( (maxOrMin != "max") && (maxOrMin != "min") ) return false;
01628 else instanceData->objectives->obj[arrayIndex]->maxOrMin = maxOrMin;
01629 instanceData->objectives->obj[arrayIndex]->constant = constant;
01630 instanceData->objectives->obj[arrayIndex]->weight = weight;
01631 int n = objectiveCoefficients->number;
01632 instanceData->objectives->obj[arrayIndex]->numberOfObjCoef = n;
01633 if(n == 0){
01634 instanceData->objectives->obj[arrayIndex]->coef = 0;
01635 }
01636 else{
01637 int i = 0;
01638 instanceData->objectives->obj[arrayIndex]->coef = new ObjCoef*[n];
01639 for(i = 0; i < n; i++) instanceData->objectives->obj[arrayIndex]->coef[i] = new ObjCoef();
01640 for(i = 0; i < n; i++){
01641 instanceData->objectives->obj[arrayIndex]->coef[i]->idx = objectiveCoefficients->indexes[i];
01642 instanceData->objectives->obj[arrayIndex]->coef[i]->value = objectiveCoefficients->values[i];
01643 }
01644 }
01645 return true;
01646 }
01647
01648 bool OSInstance::setObjectives(int number, string *names, string *maxOrMins, double *constants, double *weights, SparseVector **objectiveCoefficients){
01649 if(number < 0) return false;
01650 try{
01651 if(instanceData->objectives == NULL){
01652 throw ErrorClass("there is no objectives object");
01653 }
01654 if(instanceData->objectives->numberOfObjectives != number){
01655 throw ErrorClass("input number of objective not equal to number in class");
01656 }
01657 if(number == 0) return true;
01658 int i = 0;
01659 for(i = 0; i < number; i++)instanceData->objectives->obj[i] = new Objective();
01660 int j = 0;
01661 if(names != NULL){
01662 for(i = 0; i < number; i++) instanceData->objectives->obj[i]->name = names[i];
01663 }
01664 if(maxOrMins != NULL){
01665 for(i = 0; i < number; i++){
01666 if(maxOrMins[i] == "" || (maxOrMins[i].compare("max") != 0 && maxOrMins[i].compare("min") !=0)) return false;
01667 instanceData->objectives->obj[i]->maxOrMin = maxOrMins[i];
01668 }
01669 }
01670 if(constants != NULL){
01671 for(i = 0; i < number; i++) instanceData->objectives->obj[i]->constant = constants[i];
01672 }
01673 if(weights != NULL){
01674 for(i = 0; i < number; i++) instanceData->objectives->obj[i]->weight = weights[i];
01675 }
01676 if(objectiveCoefficients != NULL){
01677 for(i = 0; i < number; i++){
01678 int n = (&objectiveCoefficients[i] == NULL || objectiveCoefficients[i]->indexes == NULL)?0:objectiveCoefficients[i]->number;
01679 instanceData->objectives->obj[i]->numberOfObjCoef = n;
01680 if(n == 0){
01681 instanceData->objectives->obj[i]->coef = NULL;
01682 }
01683 else{
01684 instanceData->objectives->obj[i]->coef = new ObjCoef*[n];
01685 for(j = 0; j < n; j++){
01686 instanceData->objectives->obj[i]->coef[j] = new ObjCoef();
01687 instanceData->objectives->obj[i]->coef[j]->idx = objectiveCoefficients[i]->indexes[j];
01688 instanceData->objectives->obj[i]->coef[j]->value = objectiveCoefficients[i]->values[j];
01689 }
01690 }
01691 }
01692 }
01693 return true;
01694 }
01695 catch(const ErrorClass& eclass){
01696 throw ErrorClass( eclass.errormsg);
01697 }
01698 }
01699
01700
01701 bool OSInstance::setConstraintNumber(int number){
01702 if(number < 0) return false;
01703 if(instanceData->constraints == NULL) instanceData->constraints = new Constraints();
01704 if(number == 0){
01705 instanceData->constraints->numberOfConstraints = 0;
01706 instanceData->constraints->con = 0;
01707 return true;
01708 }
01709 instanceData->constraints->numberOfConstraints = number;
01710 if(instanceData->constraints->con == 0 ){
01711 instanceData->constraints->con = new Constraint*[number];
01712 }
01713 return true;
01714 }
01715
01716 bool OSInstance::addConstraint(int index, string name, double lowerBound, double upperBound, double constant) {
01717 instanceData->constraints->con[ index] = new Constraint();
01718 if(index < 0 || instanceData->constraints->numberOfConstraints <= 0 || index >= instanceData->constraints->numberOfConstraints) return false;
01719 instanceData->constraints->con[ index]->name = name;
01720 if(lowerBound != -OSDBL_MAX && lowerBound != -OSDBL_MAX) instanceData->constraints->con[ index]->lb = lowerBound;
01721 if(upperBound != OSDBL_MAX && upperBound != OSDBL_MAX)instanceData->constraints->con[ index]->ub = upperBound;
01722 instanceData->constraints->con[ index]->constant = constant;
01723 return true;
01724 }
01725
01726
01727 bool OSInstance::setConstraints(int number, string* names, double* lowerBounds, double* upperBounds, double* constants){
01728 if(number < 0) return false;
01729 if(number == 0){
01730
01731
01732
01733
01734 return true;
01735 }
01736 try{
01737
01738 if(instanceData->constraints == NULL){
01739 throw ErrorClass("there is no constraints object");
01740 }
01741 if(instanceData->constraints->numberOfConstraints != number){
01742 throw ErrorClass("input number of constrasints not equal to number in class");
01743 }
01744 int i = 0;
01745 for(i = 0; i < number; i++){
01746 instanceData->constraints->con[i] = new Constraint();
01747 }
01748 if(names != NULL){
01749 for(i = 0; i < number; i++) instanceData->constraints->con[i]->name = names[i];
01750 }
01751 if(lowerBounds != NULL){
01752 for(i = 0; i < number; i++){
01753 if(lowerBounds[i] != -OSDBL_MAX && lowerBounds[i] != -OSDBL_MAX)instanceData->constraints->con[i]->lb = lowerBounds[i];
01754 }
01755 }
01756 if(upperBounds != NULL){
01757 for(i = 0; i < number; i++){
01758 if(upperBounds[i] != OSDBL_MAX && upperBounds[i] != OSDBL_MAX)instanceData->constraints->con[i]->ub = upperBounds[i];
01759 }
01760 }
01761 if(constants != NULL){
01762 for(i = 0; i < number; i++) instanceData->constraints->con[i]->constant = constants[i];
01763 }
01764 return true;
01765 }
01766 catch(const ErrorClass& eclass){
01767 throw ErrorClass( eclass.errormsg);
01768 }
01769 }
01770
01771 bool OSInstance::setLinearConstraintCoefficients(int numberOfValues, bool isColumnMajor,
01772 double* values, int valuesBegin, int valuesEnd,
01773 int* indexes, int indexesBegin, int indexesEnd,
01774 int* starts, int startsBegin, int startsEnd){
01775 if(numberOfValues < 0) return false;
01776 if(instanceData->linearConstraintCoefficients == NULL) instanceData->linearConstraintCoefficients = new LinearConstraintCoefficients() ;
01777 if(numberOfValues == 0) return true;
01778 if((values == 0 ) ||
01779 (valuesBegin < 0 || (valuesEnd - valuesBegin + 1) != numberOfValues) ||
01780 (indexes == 0) ||
01781 (indexesBegin < 0 || (indexesEnd - indexesBegin + 1) != numberOfValues) ||
01782 (starts == 0 ) ||
01783 (startsBegin < 0 || startsBegin >= startsEnd)) return false;
01784 instanceData->linearConstraintCoefficients->numberOfValues = numberOfValues;
01785 int i = 0;
01786
01787 if(instanceData->linearConstraintCoefficients->start == NULL) instanceData->linearConstraintCoefficients->start = new IntVector();
01788 if(startsBegin == 0 ){
01789 instanceData->linearConstraintCoefficients->start->el = starts;
01790 }
01791 else{
01792 instanceData->linearConstraintCoefficients->start->el = new int[startsEnd - startsBegin + 1];
01793 int k = 0;
01794 for(i = startsBegin; i <= startsEnd; i++){
01795 instanceData->linearConstraintCoefficients->start->el[k] = starts[i];
01796 k++;
01797 }
01798 }
01799
01800 if(instanceData->linearConstraintCoefficients->value == NULL) instanceData->linearConstraintCoefficients->value = new DoubleVector();
01801 if(valuesBegin == 0 ){
01802 instanceData->linearConstraintCoefficients->value->el = values;
01803 }
01804 else{
01805 instanceData->linearConstraintCoefficients->value->el = new double[numberOfValues];
01806 int k = 0;
01807 for(i = valuesBegin; i <= valuesEnd; i++){
01808 instanceData->linearConstraintCoefficients->value->el[k] = values[i];
01809 k++;
01810 }
01811 }
01812
01813 if(instanceData->linearConstraintCoefficients->rowIdx == NULL) instanceData->linearConstraintCoefficients->rowIdx = new IntVector();
01814 if(instanceData->linearConstraintCoefficients->colIdx == NULL) instanceData->linearConstraintCoefficients->colIdx = new IntVector();
01815 if(isColumnMajor){
01816 if(indexesBegin == 0 ){
01817 instanceData->linearConstraintCoefficients->rowIdx->el = indexes;
01818 }
01819 else{
01820 instanceData->linearConstraintCoefficients->rowIdx->el = new int[numberOfValues];
01821 int k = 0;
01822 for(i = indexesBegin; i <= indexesEnd; i++){
01823 instanceData->linearConstraintCoefficients->rowIdx->el[k] = indexes[i];
01824 k++;
01825 }
01826 }
01827 }
01828 else{
01829 if(indexesBegin == 0 ){
01830 instanceData->linearConstraintCoefficients->colIdx->el = indexes;
01831 }
01832 else{
01833 instanceData->linearConstraintCoefficients->colIdx->el = new int[numberOfValues];
01834 int k = 0;
01835 for(i = indexesBegin; i <= indexesEnd; i++){
01836 instanceData->linearConstraintCoefficients->colIdx->el[k] = indexes[i];
01837 k++;
01838 }
01839 }
01840 }
01841 return true;
01842 }
01843
01844 bool OSInstance::setQuadraticTerms(int number,
01845 int* rowIndexes, int* varOneIndexes, int* varTwoIndexes, double* coefficients,
01846 int begin, int end){
01847 if(number < 0) return false;
01848 if(number != (end - begin) + 1) return false;
01849 if(number == 0){
01850 instanceData->quadraticCoefficients = 0;
01851 return true;
01852 }
01853 if( ((end - begin + 1) != number) ||
01854 (rowIndexes == 0) ||
01855 (varOneIndexes == 0) ||
01856 (varTwoIndexes == 0) ||
01857 (coefficients == 0) ) return false;
01858 instanceData->quadraticCoefficients = new QuadraticCoefficients();
01859 instanceData->quadraticCoefficients->numberOfQuadraticTerms = number;
01860 int i = 0;
01861 instanceData->quadraticCoefficients->qTerm = new QuadraticTerm*[number];
01862 for(i = 0; i < number; i++) instanceData->quadraticCoefficients->qTerm[i] = new QuadraticTerm();
01863 int k = 0;
01864 for(i = begin; i <= end; i++){
01865 instanceData->quadraticCoefficients->qTerm[k]->idx = rowIndexes[i];
01866 instanceData->quadraticCoefficients->qTerm[k]->idxOne = varOneIndexes[i];
01867 instanceData->quadraticCoefficients->qTerm[k]->idxTwo = varTwoIndexes[i];
01868 instanceData->quadraticCoefficients->qTerm[k]->coef = coefficients[i];
01869 k++;
01870 }
01871 return true;
01872 }
01873
01874 bool OSInstance::setQuadraticTermsInNonlinearExpressions(int numQPTerms, int* rowIndexes, int* varOneIndexes, int* varTwoIndexes, double* coefficients){
01875 instanceData->nonlinearExpressions->numberOfNonlinearExpressions = numQPTerms;
01876 instanceData->nonlinearExpressions->nl = new Nl*[ numQPTerms ];
01877
01878 OSnLNode *nlNodePoint;
01879 OSnLNodeVariable *nlNodeVariablePoint;
01880 std::vector<OSnLNode*> nlNodeVec;
01881
01882
01883 int i;
01884 for(i = 0; i < numQPTerms; i++){
01885 instanceData->nonlinearExpressions->nl[ i] = new Nl();
01886 instanceData->nonlinearExpressions->nl[ i]->idx = rowIndexes[ i];
01887 instanceData->nonlinearExpressions->nl[ i]->osExpressionTree = new OSExpressionTree();
01888
01889 nlNodeVariablePoint = new OSnLNodeVariable();
01890 nlNodeVariablePoint->idx = varOneIndexes[ i];
01891
01892 nlNodeVariablePoint->coef = coefficients[ i];
01893 nlNodeVec.push_back( nlNodeVariablePoint);
01894
01895 nlNodeVariablePoint = new OSnLNodeVariable();
01896 nlNodeVariablePoint->idx = varTwoIndexes[ i];
01897 nlNodeVec.push_back( nlNodeVariablePoint);
01898
01899 nlNodePoint = new OSnLNodeTimes();
01900 nlNodeVec.push_back( nlNodePoint);
01901
01902
01903 instanceData->nonlinearExpressions->nl[ i]->osExpressionTree->m_treeRoot =
01904 nlNodeVec[ 0]->createExpressionTreeFromPostfix( nlNodeVec);
01905 nlNodeVec.clear();
01906 }
01907 return true;
01908 }
01909
01910 bool OSInstance::initializeNonLinearStructures( ){
01911 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
01912 if( m_bNonLinearStructuresInitialized == true) return true;
01913 if( m_bProcessVariables == false) processVariables();
01914 if( m_bProcessObjectives == false) processObjectives();
01915 if( m_bProcessConstraints == false) processConstraints();
01916 m_iVariableNumber = instanceData->variables->numberOfVariables;
01917 m_iConstraintNumber = instanceData->constraints->numberOfConstraints;
01918 m_iObjectiveNumber = instanceData->objectives->numberOfObjectives;
01919
01920 if( m_bProcessExpressionTrees == false) getAllNonlinearExpressionTrees();
01921
01922 if( m_bDuplicateExpressionTreesMap == false) duplicateExpressionTreesMap();
01923
01924 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
01925 (posMapExpTree->second)->getVariableIndiciesMap() ;
01926 }
01927
01928 if(getNumberOfQuadraticTerms() > 0) addQTermsToExressionTree();
01929
01930 getAllNonlinearVariablesIndexMap( );
01931 getDenseObjectiveCoefficients();
01932 m_mdConstraintFunctionValues = new double[ this->instanceData->constraints->numberOfConstraints];
01933 m_mdObjectiveFunctionValues = new double[ this->instanceData->objectives->numberOfObjectives];
01934
01935 m_bNonLinearStructuresInitialized = true;
01936 return true;
01937 }
01938
01939 SparseJacobianMatrix *OSInstance::getJacobianSparsityPattern( ){
01940
01941
01942
01943
01944
01945 if( m_bSparseJacobianCalculated == true) return m_sparseJacMatrix;
01946
01947
01948 getLinearConstraintCoefficientMajor();
01949
01950 if( m_bNonLinearStructuresInitialized == false) initializeNonLinearStructures( );
01951 try{
01952 if( m_bColumnMajor == true){
01953 if( getSparseJacobianFromColumnMajor( ) == false) throw ErrorClass("An error occurred in getSpareJacobianFromColumnMajor");
01954 }
01955 else {
01956 if( getSparseJacobianFromRowMajor( ) == false) throw ErrorClass("An error occurred in getSpareJacobianFromRowMajor");
01957 }
01958 }
01959 catch(const ErrorClass& eclass){
01960 throw ErrorClass( eclass.errormsg);
01961 }
01962
01963 m_sparseJacMatrix = new SparseJacobianMatrix();
01964
01965
01966 m_sparseJacMatrix->bDeleteArrays = false;
01967 m_sparseJacMatrix->valueSize = m_iJacValueSize;
01968 m_sparseJacMatrix->starts = m_miJacStart;
01969 m_sparseJacMatrix->conVals = m_miJacNumConTerms;
01970 m_sparseJacMatrix->indexes = m_miJacIndex;
01971 m_sparseJacMatrix->values = m_mdJacValue;
01972 m_bSparseJacobianCalculated = true;
01973 return m_sparseJacMatrix;
01974 }
01975
01976 bool OSInstance::addQTermsToExressionTree(){
01977 int i, k, idx;
01978
01979 int numQTerms = instanceData->quadraticCoefficients->numberOfQuadraticTerms;
01980 if(numQTerms <= 0 || m_bQTermsAdded == true) return true;
01981 OSnLNodeVariable* nlNodeVariableOne;
01982 OSnLNodeVariable* nlNodeVariableTwo;
01983 OSnLNodeTimes* nlNodeTimes;
01984 OSnLNodePlus* nlNodePlus;
01985 OSExpressionTree* expTree;
01986 getQuadraticTerms();
01987 std::cout << "PROCESSING QUADRATIC TERMS" << std::endl;
01988 for(i = 0; i < numQTerms; i++){
01989 idx = m_quadraticTerms->rowIndexes[ i];
01990 std::cout << "PROCESSING QTERM = " << i <<std::endl;
01991
01992 if( m_mapExpressionTreesMod.find( idx) != m_mapExpressionTreesMod.end() ) {
01993
01994
01995
01996 nlNodeVariableOne = new OSnLNodeVariable();
01997 nlNodeVariableOne->idx = m_quadraticTerms->varOneIndexes[ i];
01998
01999
02000 expTree = m_mapExpressionTreesMod[ idx];
02001 if( expTree->m_bIndexMapGenerated == false) expTree->getVariableIndiciesMap();
02002 if( (*expTree->mapVarIdx).find( nlNodeVariableOne->idx) == (*expTree->mapVarIdx).end() ){
02003
02004 k = (*expTree->mapVarIdx).size();
02005 (*expTree->mapVarIdx)[ nlNodeVariableOne->idx] = k + 1;
02006 std::cout << "ADDED THE FOLLOWING VAIRABLE TO THE MAP" << nlNodeVariableOne->idx << std::endl;
02007 }
02008 nlNodeVariableOne->coef = m_quadraticTerms->coefficients[ i];
02009 nlNodeVariableTwo = new OSnLNodeVariable();
02010 nlNodeVariableTwo->idx = m_quadraticTerms->varTwoIndexes[ i];
02011
02012
02013 if( (*expTree->mapVarIdx).find( nlNodeVariableTwo->idx) == (*expTree->mapVarIdx).end() ){
02014
02015 k = (*expTree->mapVarIdx).size();
02016 (*expTree->mapVarIdx)[ nlNodeVariableTwo->idx] = k + 1;
02017 std::cout << "ADDED THE FOLLOWING VAIRABLE TO THE MAP" << nlNodeVariableTwo->idx << std::endl;
02018 }
02019 nlNodeVariableTwo->coef = 1.;
02020
02021 nlNodeTimes = new OSnLNodeTimes();
02022 nlNodeTimes->m_mChildren[ 0] = nlNodeVariableOne;
02023 nlNodeTimes->m_mChildren[ 1] = nlNodeVariableTwo;
02024
02025 nlNodePlus = new OSnLNodePlus();
02026 nlNodePlus->m_mChildren[ 0] = expTree->m_treeRoot;
02027 nlNodePlus->m_mChildren[ 1] = nlNodeTimes;
02028
02029 expTree->m_treeRoot = nlNodePlus ;
02030
02031 if(expTree->m_bIndexMapGenerated == true){
02032 delete expTree->mapVarIdx;
02033 expTree->mapVarIdx = NULL;
02034 expTree->m_bIndexMapGenerated = false;
02035 }
02036
02037
02038
02039 }
02040 else{
02041
02042 nlNodeVariableOne = new OSnLNodeVariable();
02043 nlNodeVariableOne->idx = m_quadraticTerms->varOneIndexes[ i];
02044 nlNodeVariableOne->coef = m_quadraticTerms->coefficients[ i];
02045 nlNodeVariableTwo = new OSnLNodeVariable();
02046 nlNodeVariableTwo->idx = m_quadraticTerms->varTwoIndexes[ i];
02047 nlNodeVariableTwo->coef = 1.;
02048
02049 nlNodeTimes = new OSnLNodeTimes();
02050 nlNodeTimes->m_mChildren[ 0] = nlNodeVariableOne;
02051 nlNodeTimes->m_mChildren[ 1] = nlNodeVariableTwo;
02052
02053 expTree = new OSExpressionTree();
02054 expTree->m_treeRoot = nlNodeTimes ;
02055 expTree->mapVarIdx = expTree->getVariableIndiciesMap();
02056 m_mapExpressionTreesMod[ idx ] = expTree;
02057 if(idx < 0){
02058 m_iObjectiveNumberNonlinear++;
02059 m_bProcessExpressionTrees = true;
02060 }
02061 else{
02062 m_iConstraintNumberNonlinear++;
02063 m_bProcessExpressionTrees = true;
02064 }
02065 std::cout << "NUMBER OF EXPRESSION TREES = " << m_mapExpressionTreesMod.size() <<std::endl;
02066 std::cout << "NUMBER OF NONLINEAR OBJECTIVES = " << getNumberOfNonlinearObjectives() <<std::endl;
02067 }
02068
02069 if(m_iNonlinearExpressionNumber <= 0) m_mapExpressionTrees = m_mapExpressionTreesMod;
02070 m_bQTermsAdded =true;
02071 }
02072 return true;
02073 }
02074
02075 double OSInstance::calculateFunctionValue(int idx, double *x, bool new_x){
02076 try{
02077
02078 int i, j;
02079 double dvalue = 0;
02080 if( m_binitForAlgDiff == false) initForAlgDiff();
02081 if( m_bSparseJacobianCalculated == false) getJacobianSparsityPattern();
02082 if(idx >= 0){
02083
02084 if( getConstraintNumber() <= idx ) throw
02085 ErrorClass("constraint index not valid in OSInstance::calculateFunctionValue");
02086 if( new_x == false) return *(m_mdConstraintFunctionValues + idx);
02087
02088 if( m_mapExpressionTreesMod.find( idx) != m_mapExpressionTreesMod.end() ){
02089 dvalue = m_mapExpressionTreesMod[ idx]->calculateFunction( x, new_x);
02090
02091 }
02092
02093
02094 i = m_sparseJacMatrix->starts[ idx];
02095 j = m_sparseJacMatrix->starts[ idx + 1 ];
02096
02097 while ( (i - m_sparseJacMatrix->starts[ idx]) < m_sparseJacMatrix->conVals[ idx] ){
02098
02099
02100 dvalue += m_sparseJacMatrix->values[ i]*x[ m_sparseJacMatrix->indexes[ i] ];
02101 i++;
02102 }
02103
02104 dvalue += m_mdConstraintConstants[ idx ];
02105 return dvalue;
02106 }
02107 else{
02108
02109 if( getObjectiveNumber() <= ( abs( idx) - 1) ) throw
02110 ErrorClass("objective function index not valid in OSInstance::calculateFunctionValue");
02111 if( new_x == false) return *(m_mdObjectiveFunctionValues + ( abs( idx) - 1));
02112
02113 if( m_mapExpressionTreesMod.find( idx) != m_mapExpressionTreesMod.end() ){
02114 dvalue = m_mapExpressionTreesMod[ idx]->calculateFunction( x, new_x);
02115 }
02116
02117 SparseVector **objCoef = getObjectiveCoefficients();
02118 SparseVector *obj = objCoef[ abs( idx) - 1];
02119 for(i = 0; i < obj->number; i++){
02120 dvalue += x[ obj->indexes[i]]*(obj->values[ i]);
02121 }
02122
02123 dvalue += m_mdObjectiveConstants[ abs( idx) - 1 ];
02124
02125 *(m_mdObjectiveFunctionValues + ( abs( idx) - 1)) = dvalue;
02126 return *(m_mdObjectiveFunctionValues + ( abs( idx) - 1));
02127 }
02128 }
02129 catch(const ErrorClass& eclass){
02130 throw ErrorClass( eclass.errormsg);
02131 }
02132 }
02133
02134
02135 double *OSInstance::calculateAllConstraintFunctionValues( double* x, double *objLambda, double *conLambda,
02136 bool new_x, int highestOrder){
02137 try{
02138 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated) )
02139 getIterateResults(x, objLambda, conLambda, new_x, highestOrder);
02140 }
02141 catch(const ErrorClass& eclass){
02142 throw ErrorClass( eclass.errormsg);
02143 }
02144 return m_mdConstraintFunctionValues;
02145 }
02146
02147
02148 double *OSInstance::calculateAllConstraintFunctionValues(double* x, bool new_x){
02149 try{
02150 m_iHighestOrderEvaluated = -1;
02151 if( new_x == false) return m_mdConstraintFunctionValues;
02152 int idx, numConstraints;
02153 numConstraints = getConstraintNumber();
02154
02155 for(idx = 0; idx < numConstraints; idx++){
02156 m_mdConstraintFunctionValues[ idx] = calculateFunctionValue(idx, x, new_x);
02157 }
02158
02159 }
02160 catch(const ErrorClass& eclass){
02161 throw ErrorClass( eclass.errormsg);
02162 }
02163 return m_mdConstraintFunctionValues;
02164 }
02165
02166
02167 double *OSInstance::calculateAllObjectiveFunctionValues( double* x, double *objLambda, double *conLambda,
02168 bool new_x, int highestOrder){
02169 try{
02170 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated) )
02171 getIterateResults(x, objLambda, conLambda, new_x, highestOrder);
02172 }
02173 catch(const ErrorClass& eclass){
02174 throw ErrorClass( eclass.errormsg);
02175 }
02176 return m_mdObjectiveFunctionValues;
02177 }
02178
02179
02180 double *OSInstance::calculateAllObjectiveFunctionValues( double* x, bool new_x){
02181 try{
02182 m_iHighestOrderEvaluated = -1;
02183 if( new_x == false) return m_mdObjectiveFunctionValues;
02184 int idx, numObjectives;
02185 numObjectives = getObjectiveNumber();
02186
02187 for(idx = 0; idx < numObjectives; idx++){
02188 m_mdObjectiveFunctionValues[ idx] = calculateFunctionValue(-idx -1, x, new_x);
02189 }
02190 }
02191 catch(const ErrorClass& eclass){
02192 throw ErrorClass( eclass.errormsg);
02193 }
02194 return m_mdObjectiveFunctionValues;
02195 }
02196
02197
02198 SparseJacobianMatrix *OSInstance::calculateAllConstraintFunctionGradients(double* x, double *objLambda, double *conLambda,
02199 bool new_x, int highestOrder){
02200 try{
02201 if(highestOrder < 1 ) throw ErrorClass("When calling calculateAllConstraintFunctionGradients highestOrder should be 1 or 2");
02202 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated) )
02203 getIterateResults(x, objLambda, conLambda, new_x, highestOrder);
02204 }
02205 catch(const ErrorClass& eclass){
02206 throw ErrorClass( eclass.errormsg);
02207 }
02208 return m_sparseJacMatrix;
02209 }
02210
02211
02212
02213 SparseVector *OSInstance::calculateConstraintFunctionGradient(double* x, double *objLambda, double *conLambda,
02214 int idx, bool new_x, int highestOrder){
02215 try{
02216 if(highestOrder < 1 ) throw ErrorClass("When calling calculateConstraintFunctionGradient highestOrder should be 1 or 2");
02217 if(idx > instanceData->variables->numberOfVariables )
02218 throw ErrorClass("invalid index passed to calculateConstraintFunctionGrad");
02219 SparseVector *sp;
02220 sp = new SparseVector();
02221 sp->bDeleteArrays = true;
02222 int i;
02223 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated) )
02224 getIterateResults(x, objLambda, conLambda, new_x, highestOrder);
02225 sp->number = m_miJacStart[ idx + 1] - m_miJacStart[ idx];
02226 sp->values = new double[ sp->number];
02227 sp->indexes = new int[ sp->number];
02228 for(i = 0; i < sp->number; i++){
02229 sp->values[ i] = m_mdJacValue[ m_miJacStart[ idx] + i];
02230 sp->indexes[ i] = m_miJacIndex[ m_miJacStart[ idx] + i];
02231 }
02232 return sp;
02233 }
02234 catch(const ErrorClass& eclass){
02235 throw ErrorClass( eclass.errormsg);
02236 }
02237 }
02238
02239
02240 SparseVector *OSInstance::calculateConstraintFunctionGradient(double* x, int idx, bool new_x){
02241 try{
02242 if(idx > instanceData->variables->numberOfVariables )
02243 throw ErrorClass("invalid index passed to calculateConstraintFunctionGrad");
02244 SparseVector *sp;
02245 sp = new SparseVector();
02246 sp->bDeleteArrays = true;
02247 int i;
02248 if( new_x == true || (1 > m_iHighestOrderEvaluated) )
02249 getIterateResults(x, NULL, NULL, new_x, 1);
02250 sp->number = m_miJacStart[ idx + 1] - m_miJacStart[ idx];
02251 sp->values = new double[ sp->number];
02252 sp->indexes = new int[ sp->number];
02253 for(i = 0; i < sp->number; i++){
02254 sp->values[ i] = m_mdJacValue[ m_miJacStart[ idx] + i];
02255 sp->indexes[ i] = m_miJacIndex[ m_miJacStart[ idx] + i];
02256 }
02257 return sp;
02258 }
02259 catch(const ErrorClass& eclass){
02260 throw ErrorClass( eclass.errormsg);
02261 }
02262 }
02263
02264
02265
02266 double **OSInstance::calculateAllObjectiveFunctionGradients(double* x, double *objLambda, double *conLambda,
02267 bool new_x, int highestOrder){
02268 try{
02269 if(highestOrder < 1 ) throw ErrorClass("When calling calculateAllObjectiveFunctionGradients highestOrder should be 1 or 2");
02270 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated) ) {
02271 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02272 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02273 if(posMapExpTree->first < 0){
02274 m_mmdObjGradient[ abs( posMapExpTree->first) - 1 ] = calculateObjectiveFunctionGradient(x, objLambda, conLambda,
02275 posMapExpTree->first, new_x, highestOrder);
02276 }
02277 }
02278 }
02279 }
02280 catch(const ErrorClass& eclass){
02281 throw ErrorClass( eclass.errormsg);
02282 }
02283 return m_mmdObjGradient;
02284 }
02285
02286 double *OSInstance::calculateObjectiveFunctionGradient(double* x, double *objLambda, double *conLambda,
02287 int objIdx, bool new_x, int highestOrder){
02288
02289
02290
02291 try{
02292 if(highestOrder < 1 ) throw ErrorClass("When calling calculateObjectiveFunctionGradient highestOrder should be 1 or 2");
02293 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated) ) {
02294 int domainIdx = 0;
02295 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02296 std::map<int, int>::iterator posVarIndexMap;
02297 unsigned int i;
02298 int iHighestOrderEvaluatedStore;
02299 iHighestOrderEvaluatedStore = m_iHighestOrderEvaluated;
02300 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02301
02302 if(posMapExpTree->first == objIdx){
02303 if( new_x == true ){
02304 if( m_vdX.size() > 0) m_vdX.clear();
02305 for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
02306 m_vdX.push_back( x[ posVarIndexMap->first]) ;
02307 }
02308 if( (m_bCppADFunIsCreated == false || m_bCppADMustReTape == true ) && (m_mapExpressionTreesMod.size() > 0) ) {
02309 createCppADFun( m_vdX);
02310 }
02311 }
02312
02313 if(( new_x == true ) || (m_iHighestOrderEvaluated < 0) )this->forwardAD(0, m_vdX);
02314
02315 if(( new_x == true ) || (m_iHighestOrderEvaluated < 1) ) {
02316 m_vdRangeUnitVec[ domainIdx] = 1.;
02317 m_vdYjacval = this->reverseAD(1, m_vdRangeUnitVec);
02318 for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
02319 m_mmdObjGradient[ abs( objIdx) - 1 ][ m_miNonLinearVarsReverseMap[ i]] = m_vdYjacval[ i] +
02320 m_mmdDenseObjectiveCoefficients[ abs( objIdx) - 1][ m_miNonLinearVarsReverseMap[ i]];
02321 }
02322 }
02323 m_iHighestOrderEvaluated = iHighestOrderEvaluatedStore;
02324 m_vdRangeUnitVec[ domainIdx] = 0.;
02325
02326 break;
02327 }
02328 domainIdx++;
02329 }
02330 }
02331 }
02332 catch(const ErrorClass& eclass){
02333 throw ErrorClass( eclass.errormsg);
02334 }
02335 return m_mmdObjGradient[abs( objIdx) - 1];
02336 }
02337
02338
02339 double *OSInstance::calculateObjectiveFunctionGradient(double* x, int objIdx, bool new_x){
02340 try{
02341 int domainIdx = 0;
02342 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02343 std::map<int, int>::iterator posVarIndexMap;
02344 unsigned int i;
02345 int iHighestOrderEvaluatedStore;
02346 iHighestOrderEvaluatedStore = m_iHighestOrderEvaluated;
02347 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02348 if(posMapExpTree->first == objIdx){
02349 if( new_x == true ){
02350 if( m_vdX.size() > 0) m_vdX.clear();
02351 for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
02352 m_vdX.push_back( x[ posVarIndexMap->first]) ;
02353 }
02354 if( (m_bCppADFunIsCreated == false || m_bCppADMustReTape == true ) && (m_mapExpressionTreesMod.size() > 0) ) {
02355 createCppADFun( m_vdX);
02356 }
02357 }
02358
02359 if(( new_x == true ) || (m_iHighestOrderEvaluated < 0) )this->forwardAD(0, m_vdX);
02360
02361 if(( new_x == true ) || (m_iHighestOrderEvaluated < 1) ) {
02362 m_vdRangeUnitVec[ domainIdx] = 1.;
02363 m_vdYjacval = this->reverseAD(1, m_vdRangeUnitVec);
02364 for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
02365 m_mmdObjGradient[ abs( objIdx) - 1 ][ m_miNonLinearVarsReverseMap[ i]] = m_vdYjacval[ i] +
02366 m_mmdDenseObjectiveCoefficients[ abs( objIdx) - 1][ m_miNonLinearVarsReverseMap[ i]];
02367 }
02368 }
02369 m_iHighestOrderEvaluated = iHighestOrderEvaluatedStore;
02370 m_vdRangeUnitVec[ domainIdx] = 0.;
02371
02372 break;
02373 }
02374 domainIdx++;
02375 }
02376 }
02377 catch(const ErrorClass& eclass){
02378 throw ErrorClass( eclass.errormsg);
02379 }
02380 return m_mmdObjGradient[abs( objIdx) - 1];
02381 }
02382
02383
02384 SparseHessianMatrix *OSInstance::calculateLagrangianHessian( double* x, double *objLambda, double *conLambda,
02385 bool new_x, int highestOrder){
02386 try{
02387 if(highestOrder != 2 ) throw ErrorClass("When calling calculateLagrangianHessian highestOrder should be 2");
02388 if( new_x == true || (highestOrder > m_iHighestOrderEvaluated) ) {
02389
02390 getIterateResults(x, objLambda, conLambda, new_x, highestOrder);
02391 }
02392 }
02393 catch(const ErrorClass& eclass){
02394 throw ErrorClass( eclass.errormsg);
02395 }
02396 return m_LagrangianSparseHessian;
02397 }
02398
02399 SparseHessianMatrix *OSInstance::calculateHessian(double* x, int idx, bool new_x){
02400 try{
02401 if(idx > instanceData->variables->numberOfVariables )
02402 throw ErrorClass("invalid index passed to calculateHessian");
02403 double *objLambda = new double[ getObjectiveNumber() ];
02404 double *conLambda = new double[ getConstraintNumber() ];
02405
02406
02407 int i;
02408
02409 for(i = 0; i < getObjectiveNumber(); i++){
02410 objLambda[ i] = 0.0;
02411 }
02412 for(i = 0; i < getConstraintNumber(); i++){
02413 conLambda[ i] = 0.0;
02414 }
02415
02416
02417 if(idx < 0){
02418 objLambda[ abs(idx) - 1] = 1.0;
02419 }
02420 else{
02421 conLambda[ idx] = 1.0;
02422 }
02423 if( new_x == true || (2 > m_iHighestOrderEvaluated) ) {
02424 getIterateResults(x, objLambda, conLambda, new_x, 2);
02425 std::cout << "CALL getIterateResults() FROM calculateHessian" << std::endl;
02426 }
02427 delete[] objLambda;
02428 delete[] conLambda;
02429 }
02430 catch(const ErrorClass& eclass){
02431 throw ErrorClass( eclass.errormsg);
02432 }
02433 return m_LagrangianSparseHessian;
02434 }
02435
02436
02437
02438 bool OSInstance::getSparseJacobianFromColumnMajor( ){
02439
02440 if( m_bColumnMajor == false) return false;
02441 int iNumRowStarts = getConstraintNumber() + 1;
02442 int i,j, iTemp;
02443 int iNumVariableStarts = getVariableNumber() ;
02444 int *start = NULL;
02445 int *index = NULL;
02446 double *value = NULL;
02447 if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02448 start = this->instanceData->linearConstraintCoefficients->start->el;
02449 index = this->instanceData->linearConstraintCoefficients->rowIdx->el;
02450 value = this->instanceData->linearConstraintCoefficients->value->el;
02451 }
02452 m_miJacStart = new int[ iNumRowStarts];
02453 m_miJacNumConTerms = new int[ getConstraintNumber()];
02454 OSnLNodePlus *nlNodePlus;
02455 OSnLNodeVariable *nlNodeVariable;
02456 OSExpressionTree *expTree = NULL;
02457
02458 for ( i = 0; i < iNumRowStarts; i++){
02459 m_miJacStart [ i ] = 0;
02460
02461 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
02462
02463 m_mapExpressionTreesMod[ i]->getVariableIndiciesMap();
02464
02465 }
02466 }
02467
02468 if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02469
02470 for (i = 0; i < iNumVariableStarts; i++){
02471 for (j = start[i]; j < start[ i + 1 ]; j++){
02472
02473
02474
02475
02476
02477
02478 if( (m_mapExpressionTreesMod.find( index[ j]) != m_mapExpressionTreesMod.end() ) &&
02479 ( (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).find( i) != (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).end()) ){
02480
02481
02482
02483 expTree = m_mapExpressionTreesMod[ index[j] ];
02484 nlNodeVariable = new OSnLNodeVariable();
02485 nlNodeVariable->coef = value[ j];
02486 nlNodeVariable->idx = i;
02487 nlNodePlus = new OSnLNodePlus();
02488 nlNodePlus->m_mChildren[ 0] = m_mapExpressionTreesMod[ index[ j] ]->m_treeRoot;
02489 nlNodePlus->m_mChildren[ 1] = nlNodeVariable;
02490
02491 expTree->m_treeRoot = nlNodePlus ;
02492
02493
02494
02495
02496 }
02497 else{
02498 m_miJacStart[ index[j] + 1] ++;
02499 }
02500 }
02501 }
02502 }
02503
02504
02505
02506
02507 m_miJacStart[0] = 0;
02508 for (i = 1; i < iNumRowStarts; i++ ){
02509 m_miJacNumConTerms[ i - 1] = m_miJacStart[i];
02510 if( m_mapExpressionTreesMod.find( i - 1) != m_mapExpressionTreesMod.end() ){
02511 m_miJacStart[i] += (m_miJacStart[i - 1] + (*m_mapExpressionTreesMod[ i - 1]->mapVarIdx).size() );
02512 }
02513 else{
02514 m_miJacStart[i] += m_miJacStart[i - 1];
02515 }
02516 }
02517
02518 m_iJacValueSize = m_miJacStart[ iNumRowStarts - 1];
02519 m_miJacIndex = new int[ m_iJacValueSize];
02520 m_mdJacValue = new double[ m_iJacValueSize ];
02521
02522 if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02523
02524 for (i = 0; i < iNumVariableStarts; i++){
02525
02526
02527 for (j = start[i]; j < start[ i + 1 ]; j++){
02528
02529
02530 if( (m_mapExpressionTreesMod.find( index[ j]) == m_mapExpressionTreesMod.end() ) ||
02531 ( (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).find( i) == (*m_mapExpressionTreesMod[ index[ j]]->mapVarIdx).end()) ){
02532 iTemp = m_miJacStart[ index[j]];
02533 m_miJacIndex[ iTemp] = i;
02534 m_mdJacValue[ iTemp] = value[j];
02535 m_miJacStart[ index[j]]++;
02536 }
02537 }
02538 }
02539 }
02540
02541 std::map<int, int>::iterator posVarIdx;
02542
02543 for (i = 0; i < iNumRowStarts - 1; i++ ){
02544 m_miJacStart[ i] = m_miJacStart [ i] - m_miJacNumConTerms[ i] ;
02545 iTemp = m_miJacStart[ i] + m_miJacNumConTerms[ i];
02546
02547 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
02548 for(posVarIdx = (*m_mapExpressionTreesMod[ i]->mapVarIdx).begin(); posVarIdx
02549 != (*m_mapExpressionTreesMod[ i]->mapVarIdx).end(); ++posVarIdx){
02550 m_miJacIndex[ iTemp] = posVarIdx->first;
02551 m_mdJacValue[ iTemp] = 0;
02552 iTemp++;
02553 }
02554 }
02555 }
02556 #ifdef DEBUG
02557 std::cout << "HERE ARE ROW STARTS:" << std::endl;
02558 for (i = 0; i < iNumRowStarts; i++ ){
02559 std::cout << m_miJacStart[ i] << " ";
02560 }
02561 std::cout << std::endl << std::endl;
02562 std::cout << "HERE ARE VARIABLE INDICIES:" << std::endl;
02563 for (i = 0; i < m_miJacStart[ iNumRowStarts - 1]; i++ ){
02564 std::cout << m_miJacIndex[ i] << " ";
02565 }
02566 std::cout << std::endl << std::endl;
02567 std::cout << "HERE ARE VALUES:" << std::endl;
02568 for (i = 0; i < m_miJacStart[ iNumRowStarts - 1]; i++ ){
02569 std::cout << m_mdJacValue[ i] << " ";
02570 }
02571 std::cout << std::endl << std::endl;
02572
02573 std::cout << "HERE ARE NUMBER OF CONSTANT TERMS:" << std::endl;
02574 for (i = 0; i < iNumRowStarts - 1; i++ ){
02575 std::cout << m_miJacNumConTerms[ i ] << " ";
02576 }
02577 std::cout << std::endl << std::endl;
02578 #endif
02579 return true;
02580 }
02581
02582
02583
02584 bool OSInstance::getSparseJacobianFromRowMajor( ){
02585
02586 if( m_bColumnMajor == true) return false;
02587 int iNumJacRowStarts = getConstraintNumber() + 1;
02588 std::map<int, int>::iterator posVarIdx;
02589 int i,j, k;
02590 int *start = NULL;
02591 int *index = NULL;
02592 double *value = NULL;
02593 if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02594 start = this->instanceData->linearConstraintCoefficients->start->el;
02595 index = this->instanceData->linearConstraintCoefficients->colIdx->el;
02596 value = this->instanceData->linearConstraintCoefficients->value->el;
02597 }
02598 m_miJacStart = new int[ iNumJacRowStarts];
02599 m_miJacNumConTerms = new int[ getConstraintNumber()];
02600 OSnLNodePlus *nlNodePlus;
02601 OSnLNodeVariable *nlNodeVariable;
02602
02603
02604 for ( i = 0; i < iNumJacRowStarts; i++){
02605 m_miJacStart [ i ] = 0;
02606
02607 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
02608
02609 m_mapExpressionTreesMod[ i]->getVariableIndiciesMap();
02610
02611 }
02612 }
02613 int loopLimit = getConstraintNumber();
02614
02615
02616 if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02617
02618 for (i = 0; i < loopLimit; i++){
02619 m_miJacNumConTerms[ i] = 0;
02620 for (j = start[i]; j < start[ i + 1 ]; j++){
02621
02622
02623
02624 if( (m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ) &&
02625 ( (*m_mapExpressionTreesMod[ i]->mapVarIdx).find( index[ j]) != (*m_mapExpressionTreesMod[ i]->mapVarIdx).end()) ){
02626
02627
02628
02629 nlNodeVariable = new OSnLNodeVariable();
02630 nlNodeVariable->coef = value[ j];
02631 nlNodeVariable->idx = index[ j];
02632 nlNodePlus = new OSnLNodePlus();
02633 nlNodePlus->m_mChildren[ 0] = m_mapExpressionTreesMod[ i ]->m_treeRoot;
02634 nlNodePlus->m_mChildren[ 1] = nlNodeVariable;
02635
02636
02637
02638
02639 m_mapExpressionTreesMod[ i ]->m_treeRoot = nlNodePlus;
02640 }
02641 else{
02642
02643 m_miJacNumConTerms[ i]++;
02644 }
02645 }
02646 }
02647 }
02648
02649 m_miJacStart[0] = 0;
02650 for (i = 1; i < iNumJacRowStarts; i++ ){
02651 if( m_mapExpressionTreesMod.find( i - 1) != m_mapExpressionTreesMod.end() ){
02652 m_miJacStart[i] = m_miJacStart[i - 1] + (m_miJacNumConTerms[ i - 1] + (*m_mapExpressionTreesMod[ i - 1]->mapVarIdx).size() );
02653 }
02654 else{
02655 m_miJacStart[i] = m_miJacStart[i - 1] + m_miJacNumConTerms[ i - 1];
02656 }
02657 }
02658
02659
02660 m_iJacValueSize = m_miJacStart[ iNumJacRowStarts - 1];
02661 m_miJacIndex = new int[ m_iJacValueSize];
02662 m_mdJacValue = new double[ m_iJacValueSize ];
02663
02664
02665 if(this->instanceData->linearConstraintCoefficients->numberOfValues > 0){
02666 for (i = 0; i < loopLimit; i++){
02667 k = 0;
02668 for (j = start[i]; j < start[ i + 1 ]; j++){
02669 if( (m_mapExpressionTreesMod.find( i) == m_mapExpressionTreesMod.end() ) ||
02670 ( (*m_mapExpressionTreesMod[ i]->mapVarIdx).find( index[ j]) == (*m_mapExpressionTreesMod[ i]->mapVarIdx).end()) ){
02671 m_miJacIndex[ m_miJacStart[i] + k ] = index[ j];
02672 m_mdJacValue[ m_miJacStart[i] + k ] = value[ j];
02673 k++;
02674 }
02675 }
02676 }
02677 }
02678
02679 for (i = 0; i < loopLimit; i++ ){
02680 k = m_miJacStart[i] + m_miJacNumConTerms[ i ];
02681
02682 if( m_mapExpressionTreesMod.find( i) != m_mapExpressionTreesMod.end() ){
02683 for(posVarIdx = (*m_mapExpressionTreesMod[ i]->mapVarIdx).begin(); posVarIdx
02684 != (*m_mapExpressionTreesMod[ i]->mapVarIdx).end(); ++posVarIdx){
02685 m_miJacIndex[ k] = posVarIdx->first;
02686 m_mdJacValue[ k] = 0;
02687 k++;
02688 }
02689 }
02690 }
02691 #ifdef DEBUG
02692 std::cout << "HERE ARE ROW STARTS:" << std::endl;
02693 for (i = 0; i < iNumJacRowStarts; i++ ){
02694 std::cout << m_miJacStart[ i] << " ";
02695 }
02696 std::cout << std::endl << std::endl;
02697 std::cout << "HERE ARE VARIABLE INDICIES:" << std::endl;
02698 for (i = 0; i < m_miJacStart[ iNumJacRowStarts - 1]; i++ ){
02699 std::cout << m_miJacIndex[ i] << " ";
02700 }
02701 std::cout << std::endl << std::endl;
02702 std::cout << "HERE ARE VALUES:" << std::endl;
02703 for (i = 0; i < m_miJacStart[ iNumJacRowStarts - 1]; i++ ){
02704 std::cout << m_mdJacValue[ i] << " ";
02705 }
02706 std::cout << std::endl << std::endl;
02707
02708 std::cout << "HERE ARE NUMBER OF CONSTANT TERMS:" << std::endl;
02709 for (i = 0; i < iNumJacRowStarts - 1; i++ ){
02710 std::cout << m_miJacNumConTerms[ i ] << " ";
02711 }
02712 std::cout << std::endl << std::endl;
02713 #endif
02714 return true;
02715 }
02716
02717 OSExpressionTree* OSInstance::getLagrangianExpTree( ){
02718 if( m_bLagrangianExpTreeCreated == true) return m_LagrangianExpTree;
02719
02720
02721
02722 if(m_bSparseJacobianCalculated == false) getJacobianSparsityPattern( );
02723 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02724 OSnLNodeTimes* nlNodeTimes = NULL;
02725 OSnLNodeVariable* nlNodeVariable = NULL;
02726 OSnLNodeSum* nlNodeSum = NULL;
02727 int numChildren = 0;
02728 int rowIdx;
02729
02730 nlNodeSum = new OSnLNodeSum();
02731 nlNodeSum->inumberOfChildren = m_mapExpressionTreesMod.size();
02732 nlNodeSum->m_mChildren = new OSnLNode*[ nlNodeSum->inumberOfChildren];
02733
02734 m_LagrangianExpTree = new OSExpressionTree();
02735 m_LagrangianExpTree->m_treeRoot = nlNodeSum;
02736
02737 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02738
02739 nlNodeVariable = new OSnLNodeVariable();
02740 nlNodeVariable->coef = 1.;
02741
02742
02743 rowIdx = posMapExpTree->first;
02744 if(rowIdx >= 0){
02745 nlNodeVariable->idx = instanceData->variables->numberOfVariables + rowIdx;
02746 }
02747 else{
02748 nlNodeVariable->idx = instanceData->variables->numberOfVariables +
02749 instanceData->constraints->numberOfConstraints + (abs(rowIdx) - 1);
02750 }
02751
02752 nlNodeTimes = new OSnLNodeTimes();
02753 nlNodeTimes->m_mChildren[ 0] = nlNodeVariable;
02754 nlNodeTimes->m_mChildren[ 1] = m_mapExpressionTreesMod[ posMapExpTree->first ]->m_treeRoot;
02755
02756 nlNodeSum->m_mChildren[ numChildren] = nlNodeTimes;
02757 numChildren++;
02758 }
02759
02760 m_LagrangianExpTree->getVariableIndiciesMap();
02761
02762
02763
02764
02765 m_bLagrangianExpTreeCreated = true;
02766 return m_LagrangianExpTree;
02767 }
02768
02769 std::map<int, int> OSInstance::getAllNonlinearVariablesIndexMap( ){
02770 if(m_bAllNonlinearVariablesIndex == true) return m_mapAllNonlinearVariablesIndex;
02771
02772
02773 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02774 std::map<int, int>::iterator posVarIdx;
02775 OSExpressionTree *expTree;
02776 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02777
02778 expTree = posMapExpTree->second;
02779 if(expTree->m_bIndexMapGenerated == false)expTree->getVariableIndiciesMap();
02780 for(posVarIdx = (*expTree->mapVarIdx).begin(); posVarIdx != (*expTree->mapVarIdx).end(); ++posVarIdx){
02781 if( m_mapAllNonlinearVariablesIndex.find( posVarIdx->first) == m_mapAllNonlinearVariablesIndex.end() ){
02782
02783 m_mapAllNonlinearVariablesIndex[ posVarIdx->first] = 1;
02784 }
02785 }
02786 }
02787 m_miNonLinearVarsReverseMap = new int[m_mapAllNonlinearVariablesIndex.size()];
02788
02789 int kount = 0;
02790
02791 for(posVarIdx = m_mapAllNonlinearVariablesIndex.begin(); posVarIdx !=m_mapAllNonlinearVariablesIndex.end(); ++posVarIdx){
02792 posVarIdx->second = kount;
02793 m_miNonLinearVarsReverseMap[ kount] = posVarIdx->first;
02794 kount++;
02795 #ifdef DEBUG
02796 std::cout << "POSITION FIRST = " << posVarIdx->first ;
02797 std::cout << " POSITION SECOND = " << posVarIdx->second << std::endl;
02798 #endif
02799 }
02800 m_iNumberOfNonlinearVariables = kount;
02801
02802 m_bAllNonlinearVariablesIndex = true;
02803 return m_mapAllNonlinearVariablesIndex;
02804 }
02805
02806 SparseHessianMatrix* OSInstance::getLagrangianHessianSparsityPattern( ){
02807
02808 if( m_bLagrangianSparseHessianCreated == true) return m_LagrangianSparseHessian;
02809 if( m_iNumberOfNonlinearVariables == 0) return NULL;
02810 if( m_binitForAlgDiff == false ) initForAlgDiff();
02811 unsigned int i = 0;
02812 int numNonz = 0;
02813
02814
02815 std::vector<double> vx;
02816 std::map<int, int>::iterator posMap1, posMap2;
02817 if( (m_bCppADFunIsCreated == false || m_bCppADMustReTape == true ) && (m_mapExpressionTreesMod.size() > 0) ) {
02818 for(posMap1 = m_mapAllNonlinearVariablesIndex.begin(); posMap1 != m_mapAllNonlinearVariablesIndex.end(); ++posMap1){
02819 vx.push_back( 1.0) ;
02820 }
02821 createCppADFun( vx);
02822 }
02823
02824
02825 std::vector<bool> r(m_iNumberOfNonlinearVariables * m_iNumberOfNonlinearVariables);
02826 unsigned int j;
02827 for(i = 0; i < m_iNumberOfNonlinearVariables; i++) {
02828 for(j = 0; j < m_iNumberOfNonlinearVariables; j++)
02829 r[ i * m_iNumberOfNonlinearVariables + j ] = false;
02830 r[ i * m_iNumberOfNonlinearVariables + i] = true;
02831 }
02832
02833 (*Fad).ForSparseJac(m_iNumberOfNonlinearVariables, r);
02834
02835
02836 unsigned int m = m_mapExpressionTreesMod.size();
02837 std::vector<bool> e( m);
02838
02839 for(i = 0; i < m; i++) e[i] = true;
02840 std::cout << "Computing Sparse Hessian" << std::endl;
02841
02842 m_vbLagHessNonz = (*Fad).RevSparseHes(m_iNumberOfNonlinearVariables, e);
02843 for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
02844
02845 for(j = i; j < m_iNumberOfNonlinearVariables; j++){
02846 if(m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j] == true) numNonz++;
02847
02848 }
02849
02850 }
02851
02852 i = 0;
02853
02854 m_LagrangianSparseHessian = new SparseHessianMatrix();
02855 m_LagrangianSparseHessian->bDeleteArrays = true;
02856 m_LagrangianSparseHessian->hessDimension = numNonz;
02857 m_LagrangianSparseHessian->hessRowIdx = new int[m_LagrangianSparseHessian->hessDimension];
02858 m_LagrangianSparseHessian->hessColIdx = new int[m_LagrangianSparseHessian->hessDimension];
02859 m_LagrangianSparseHessian->hessValues = new double[m_LagrangianSparseHessian->hessDimension];
02860
02861 numNonz = 0;
02862 for(posMap1 = m_mapAllNonlinearVariablesIndex.begin(); posMap1 != m_mapAllNonlinearVariablesIndex.end(); ++posMap1){
02863
02864 j = i;
02865 for(posMap2 = posMap1; posMap2 != m_mapAllNonlinearVariablesIndex.end(); ++posMap2){
02866
02867 if(m_vbLagHessNonz[ i*m_iNumberOfNonlinearVariables + j] == true){
02868 *(m_LagrangianSparseHessian->hessRowIdx + numNonz) = posMap1->first;
02869 *(m_LagrangianSparseHessian->hessColIdx + numNonz) = posMap2->first;
02870 numNonz++;
02871 }
02872
02873 j++;
02874 }
02875 i++;
02876 }
02877 #ifdef DEBUG
02878 std::cout << "HESSIAN SPARSITY PATTERN" << std::endl;
02879 int kj;
02880 for(kj = 0; kj < m_LagrangianSparseHessian->hessDimension; kj++){
02881 std::cout << "Row Index = " << *(m_LagrangianSparseHessian->hessRowIdx + kj) << std::endl;
02882 std::cout << "Column Index = " << *(m_LagrangianSparseHessian->hessColIdx + kj) << std::endl;
02883 }
02884 #endif
02885
02886 m_bLagrangianSparseHessianCreated = true;
02887 return m_LagrangianSparseHessian;
02888 }
02889
02890
02891 void OSInstance::duplicateExpressionTreesMap(){
02892
02893
02894 if(m_bDuplicateExpressionTreesMap == false){
02895
02896 if( m_bProcessExpressionTrees == false) getAllNonlinearExpressionTrees();
02897 m_mapExpressionTreesMod = m_mapExpressionTrees;
02898 m_bDuplicateExpressionTreesMap = true;
02899 return;
02900 }
02901 else{
02902 return;
02903 }
02904 }
02905
02906
02907 bool OSInstance::createCppADFun(std::vector<double> vdX){
02908 if(m_bCppADFunIsCreated == true) return true;
02909
02910 if(m_binitForAlgDiff == false) initForAlgDiff();
02911
02912
02913 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
02914 unsigned int i;
02915 size_t n = vdX.size();
02916
02917 CppADvector< AD<double> > vdaX( n );
02918 for(i = 0; i < n; i++){
02919 vdaX[ i] = vdX[ i];
02920
02921 }
02922
02923 CppAD::Independent( vdaX);
02929 CppAD::vector< AD<double> > m_vFG;
02930 int kount = 0;
02931 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
02932 m_vFG.push_back( (posMapExpTree->second)->m_treeRoot->constructCppADTape(&m_mapAllNonlinearVariablesIndex, &vdaX) );
02933
02934 if( m_mapCppADFunRangeIndex.find( posMapExpTree->first) == m_mapCppADFunRangeIndex.end() ){
02935
02936 m_mapCppADFunRangeIndex[ posMapExpTree->first] = kount;
02937 kount++;
02938 }
02939 }
02940
02941 std::cout << "create the function and stop recording" << std::endl;
02942 Fad = new CppAD::ADFun<double>(vdaX, m_vFG);
02943 std::cout << "range space dimension = " << m_vFG.size() << std::endl;
02944
02945 m_iHighestTaylorCoeffOrder = -1;
02946 m_bCppADFunIsCreated = true;
02947 return true;
02948 }
02949
02950
02951 std::vector<double> OSInstance::forwardAD(int p, std::vector<double> vdX){
02952 try{
02953
02954 if(m_bCppADFunIsCreated == false) createCppADFun( vdX);
02955 if(p > (m_iHighestTaylorCoeffOrder + 1) ) throw
02956 ErrorClass( "trying to calculate a p order forward when p-1 Taylor coefficient not available");
02957
02958 m_iHighestTaylorCoeffOrder = p;
02959 m_iHighestOrderEvaluated = p;
02960
02961
02962
02963 return (*Fad).Forward(p, vdX);
02964 }
02965 catch(const ErrorClass& eclass){
02966 throw ErrorClass( eclass.errormsg);
02967 }
02968 }
02969
02970
02971 std::vector<double> OSInstance::reverseAD(int p, std::vector<double> vdlambda){
02972 try{
02973 if(p == 0) throw
02974 ErrorClass( "reverseAD must have p >= 1");
02975 if(p > (m_iHighestTaylorCoeffOrder + 1) ) throw
02976 ErrorClass( "trying to calculate a p order reverse when p-1 Taylor coefficient not available");
02977
02978
02979
02980 m_iHighestOrderEvaluated = p;
02981 return (*Fad).Reverse(p, vdlambda);
02982 }
02983 catch(const ErrorClass& eclass){
02984 throw ErrorClass( eclass.errormsg);
02985 }
02986 }
02987
02988 bool OSInstance::getIterateResults( double *x, double *objLambda, double* conMultipliers,
02989 bool new_x, int highestOrder){
02990 try{
02991 if( m_binitForAlgDiff == false) initForAlgDiff();
02992 std::map<int, int>::iterator posVarIndexMap;
02993
02994 if(new_x == true){
02995 if( m_vdX.size() > 0) m_vdX.clear();
02996 for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
02997 m_vdX.push_back( x[ posVarIndexMap->first]) ;
02998 }
02999 if( (m_bCppADFunIsCreated == false || m_bCppADMustReTape == true ) && (m_mapExpressionTreesMod.size() > 0) ) {
03000 createCppADFun( m_vdX);
03001 }
03002 }
03003 else{
03004 if( m_vdX.size() == 0) {
03005 for(posVarIndexMap = m_mapAllNonlinearVariablesIndex.begin(); posVarIndexMap != m_mapAllNonlinearVariablesIndex.end(); ++posVarIndexMap){
03006 m_vdX.push_back( x[ posVarIndexMap->first]) ;
03007 }
03008 if( (m_bCppADFunIsCreated == false || m_bCppADMustReTape == true ) && (m_mapExpressionTreesMod.size() > 0) ) {
03009 createCppADFun( m_vdX);
03010 }
03011 }
03012 }
03013 switch( highestOrder){
03014 case 0:
03015 if(new_x == true || m_iHighestOrderEvaluated < 0){
03016 if(bUseExpTreeForFunEval == true){
03017 calculateAllConstraintFunctionValues( x, new_x);
03018 calculateAllObjectiveFunctionValues( x, new_x);
03019 }
03020 else{
03021 getZeroOrderResults(x, objLambda, conMultipliers);
03022 }
03023
03024 }
03025 break;
03026 case 1:
03027 if(new_x == true || m_iHighestOrderEvaluated < 0)
03028 getZeroOrderResults(x, objLambda, conMultipliers);
03029 if(new_x == true || m_iHighestOrderEvaluated < 1)
03030 getFirstOrderResults(x, objLambda, conMultipliers);
03031 break;
03032 case 2:
03033 if(new_x == true || m_iHighestOrderEvaluated < 0)
03034 getZeroOrderResults(x, objLambda, conMultipliers);
03035 if(new_x == true || m_iHighestOrderEvaluated < 2)
03036 getSecondOrderResults(x, objLambda, conMultipliers);
03037 break;
03038 default:
03039 throw ErrorClass("Derivative should be order 0, 1, or 2");
03040 }
03041 return true;
03042 }
03043 catch(const ErrorClass& eclass){
03044 throw ErrorClass( eclass.errormsg);
03045 }
03046 }
03047
03048
03049 bool OSInstance::getZeroOrderResults(double *x, double *objLambda, double *conMultipliers){
03050 try{
03051
03052 int i, j, rowNum, objNum;
03053 if( m_mapExpressionTreesMod.size() > 0){
03054 m_vdYval = this->forwardAD(0, m_vdX);
03055 }
03056
03057 for(rowNum = 0; rowNum < m_iConstraintNumber; rowNum++){
03058 m_mdConstraintFunctionValues[ rowNum] = 0.0;
03059 if( m_mapExpressionTreesMod.find( rowNum) != m_mapExpressionTreesMod.end() ){
03060 m_mdConstraintFunctionValues[ rowNum] = m_vdYval[ m_mapCppADFunRangeIndex[ rowNum]];
03061 }
03062
03063
03064 i = m_sparseJacMatrix->starts[ rowNum];
03065 j = m_sparseJacMatrix->starts[ rowNum + 1 ];
03066 while ( (i - m_sparseJacMatrix->starts[ rowNum]) < m_sparseJacMatrix->conVals[ rowNum] ){
03067 m_mdConstraintFunctionValues[ rowNum] += m_sparseJacMatrix->values[ i]*x[ m_sparseJacMatrix->indexes[ i] ];
03068 i++;
03069 }
03070
03071 m_mdConstraintFunctionValues[ rowNum] += m_mdConstraintConstants[ rowNum ];
03072 #ifdef DEBUG
03073 std::cout << "Constraint " << rowNum << " function value = " << m_mdConstraintFunctionValues[ rowNum ] << std::endl;
03074 #endif
03075 }
03076
03077 for(objNum = 0; objNum < m_iObjectiveNumber; objNum++){
03078 m_mdObjectiveFunctionValues[ objNum] = 0.0;
03079 if( m_mapExpressionTreesMod.find( -objNum -1) != m_mapExpressionTreesMod.end() ){
03080 m_mdObjectiveFunctionValues[ objNum] = m_vdYval[ objNum];
03081 }
03082 for(i = 0; i < m_iVariableNumber; i++){
03083 m_mdObjectiveFunctionValues[ objNum] += m_mmdDenseObjectiveCoefficients[ objNum][i]*x[ i];
03084 }
03085 #ifdef DEBUG
03086 std::cout << "Objective " << objNum << " function value = " << m_mdObjectiveFunctionValues[ objNum] << std::endl;
03087 #endif
03088 }
03089 return true;
03090 }
03091 catch(const ErrorClass& eclass){
03092 throw ErrorClass( eclass.errormsg);
03093 }
03094 }
03095
03096
03097
03098 bool OSInstance::getFirstOrderResults(double *x, double *objLambda, double *conMultipliers){
03099 try{
03100
03101 unsigned int i, j;
03102 int rowNum, jacIndex;
03103 unsigned int jstart, jend;
03104 int idx;
03105 OSExpressionTree *expTree = NULL;
03106 int domainIdx = 0;
03107 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03108 std::map<int, int>::iterator posVarIdx;
03109
03115 if(m_iNumberOfNonlinearVariables >= m_mapExpressionTreesMod.size() ){
03116
03117
03118 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03119 idx = posMapExpTree->first;
03120
03121 if(idx >= 0){
03122 m_vdRangeUnitVec[ domainIdx] = 1.;
03123 m_mapExpressionTreesMod[ idx]->getVariableIndiciesMap();
03124 m_vdYjacval = this->reverseAD(1, m_vdRangeUnitVec);
03125
03126 jstart = m_miJacStart[ idx] + m_miJacNumConTerms[ idx];
03127 jend = m_miJacStart[ idx + 1 ];
03128 if( (*m_mapExpressionTreesMod[ idx]->mapVarIdx).size() != (jend - jstart)) throw
03129 ErrorClass("number of partials not consistent");
03130 j = 0;
03131 jacIndex = 0;
03132 for(posVarIdx = m_mapAllNonlinearVariablesIndex.begin(); posVarIdx
03133 != m_mapAllNonlinearVariablesIndex.end(); ++posVarIdx){
03134
03135
03136
03137
03138
03139 if( (*m_mapExpressionTreesMod[ idx]->mapVarIdx).find( posVarIdx->first) != (*m_mapExpressionTreesMod[ idx]->mapVarIdx).end()){
03140 m_mdJacValue[ jstart] = m_vdYjacval[ jacIndex];
03141 jstart++;
03142 j++;
03143 }
03144 jacIndex++;
03145 }
03146
03147 m_vdRangeUnitVec[ domainIdx] = 0.;
03148 domainIdx++;
03149 }
03150 else{
03151
03152 domainIdx++;
03153 }
03154 }
03155 }
03156 else{
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
03179 m_mdJacValue[ jstart + jacIndex] = m_vdYjacval[m_iObjectiveNumberNonlinear + rowNum];
03180 }
03181 rowNum++;
03182 }
03183
03184 }
03185
03186 m_vdDomainUnitVec[i] = 0.;
03187 }
03188 }
03189 #ifdef DEBUG
03190 int k;
03191 std::cout << "JACOBIAN DATA " << std::endl;
03192 for(idx = 0; idx < m_iConstraintNumber; idx++){
03193 for(k = *(m_sparseJacMatrix->starts + idx); k < *(m_sparseJacMatrix->starts + idx + 1); k++){
03194 std::cout << "row idx = " << idx << " col idx = "<< *(m_sparseJacMatrix->indexes + k)
03195 << " value = " << *(m_sparseJacMatrix->values + k) << std::endl;
03196 }
03197 }
03198
03199 #endif
03200 return true;
03201 }
03202 catch(const ErrorClass& eclass){
03203 throw ErrorClass( eclass.errormsg);
03204 }
03205 }
03206
03207
03208 bool OSInstance::getSecondOrderResults(double *x, double *objLambda, double *conMultipliers){
03209 try{
03210
03211 unsigned int i, j;
03212 int rowNum, jacIndex;
03213 int jstart, idx;
03214 OSExpressionTree *expTree = NULL;
03215 int hessValuesIdx = 0;
03216 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03217 std::map<int, int>::iterator posVarIndexMap;
03218 if( conMultipliers == NULL) throw ErrorClass("cannot have a null vector of lagrange multipliers when calculating Hessian of Lagrangian");
03219 if( m_vdLambda.size() > 0) m_vdLambda.clear();
03220 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03221 if( posMapExpTree->first >= 0){
03222 m_vdLambda.push_back( conMultipliers[ posMapExpTree->first]);
03223 }
03224 else{
03225
03226 m_vdLambda.push_back( objLambda[ abs(posMapExpTree->first) - 1] );
03227 }
03228 }
03229 for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
03230 m_vdDomainUnitVec[i] = 1.;
03231 rowNum = 0;
03232 if( m_mapExpressionTreesMod.size() > 0){
03233 m_vdYjacval = this->forwardAD(1, m_vdDomainUnitVec);
03234 }
03235
03236
03237 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03238 idx = posMapExpTree->first;
03239
03240 if(idx >= 0){
03241
03242
03243
03244
03245 expTree = m_mapExpressionTreesMod[ idx];
03246 if( (*expTree->mapVarIdx).find( m_miNonLinearVarsReverseMap[ i]) != (*expTree->mapVarIdx).end() ){
03247 jacIndex = (*m_mapExpressionTreesMod[ idx]->mapVarIdx)[ m_miNonLinearVarsReverseMap[ i]];
03248 jstart = m_miJacStart[ idx] + m_miJacNumConTerms[ idx];
03249 m_mdJacValue[ jstart + jacIndex] = m_vdYjacval[m_iObjectiveNumberNonlinear + rowNum];
03250 }
03251 rowNum++;
03252 }
03253 else{
03254
03255
03256 m_mmdObjGradient[ (abs( idx) - 1)][ m_miNonLinearVarsReverseMap[ i]] = m_vdYjacval[ (abs( idx) - 1)] +
03257 m_mmdDenseObjectiveCoefficients[ (abs( idx) - 1)][ m_miNonLinearVarsReverseMap[ i]];
03258 }
03259 }
03260
03261 if( m_mapExpressionTreesMod.size() > 0){
03262 m_vdw = reverseAD(2, m_vdLambda);
03263 }
03264 for(j = i; j < m_iNumberOfNonlinearVariables; j++){
03265 if( m_vbLagHessNonz[i*m_iNumberOfNonlinearVariables + j] == true){
03266 m_LagrangianSparseHessian->hessValues[ hessValuesIdx] = m_vdw[ j*2 + 1];
03267 #ifdef DEBUG
03268 std::cout << "reverse 2 " << m_LagrangianSparseHessian->hessValues[ hessValuesIdx] << std::endl;
03269 #endif
03270 hessValuesIdx++;
03271 }
03272 }
03273
03274
03275 m_vdDomainUnitVec[i] = 0.;
03276 }
03277
03278 #ifdef DEBUG
03279 int k;
03280 std::cout << "JACOBIAN DATA " << std::endl;
03281 for(idx = 0; idx < m_iConstraintNumber; idx++){
03282 for(k = *(m_sparseJacMatrix->starts + idx); k < *(m_sparseJacMatrix->starts + idx + 1); k++){
03283 std::cout << "row idx = " << idx << " col idx = "<< *(m_sparseJacMatrix->indexes + k)
03284 << " value = " << *(m_sparseJacMatrix->values + k) << std::endl;
03285 }
03286 }
03287 #endif
03288 return true;
03289 }
03290 catch(const ErrorClass& eclass){
03291 throw ErrorClass( eclass.errormsg);
03292 }
03293 }
03294
03295 bool OSInstance::initForAlgDiff(){
03296 if( m_binitForAlgDiff == true ) return true;
03297 initializeNonLinearStructures( );
03298 initObjGradients();
03299 getAllNonlinearVariablesIndexMap( );
03300 if(m_bSparseJacobianCalculated == false) getJacobianSparsityPattern();
03301
03302
03303 std::map<int, OSExpressionTree*>::iterator posMapExpTree;
03304 for(posMapExpTree = m_mapExpressionTreesMod.begin(); posMapExpTree != m_mapExpressionTreesMod.end(); ++posMapExpTree){
03305 if(posMapExpTree->second->bCppADMustReTape == true) m_bCppADMustReTape = true;
03306 }
03307
03308 #ifdef DEBUG
03309 std::cout << "RETAPE == " << m_bCppADMustReTape << std::endl;
03310 #endif
03311 unsigned int i;
03312 for(i = 0; i < m_iNumberOfNonlinearVariables; i++){
03313 m_vdDomainUnitVec.push_back( 0.0 );
03314 }
03315 for(i = 0; i < m_mapExpressionTreesMod.size(); i++){
03316 m_vdRangeUnitVec.push_back( 0.0 );
03317 }
03318 m_binitForAlgDiff = true;
03319 return true;
03320 }
03321
03322 bool OSInstance::initObjGradients(){
03323 int i, j;
03324 int m, n;
03325 m = getObjectiveNumber();
03326 n = getVariableNumber();
03327 getDenseObjectiveCoefficients();
03328 m_mmdObjGradient = new double*[m];
03329 for(i = 0; i < m; i++){
03330 m_mmdObjGradient[i] = new double[n];
03331 for(j = 0; j < n; j++){
03332 m_mmdObjGradient[i][j] = m_mmdDenseObjectiveCoefficients[ i][j];
03333 #ifdef DEBUG
03334 std::cout << "m_mmdObjGradient[i][j] = " << m_mmdObjGradient[i][j] << std::endl;
03335 #endif
03336 }
03337 }
03338 return true;
03339 }