00001
00020 #include <iomanip>
00021 #include "SmiScnModel.hpp"
00022 #include "flopc.hpp"
00023
00024
00025 #include "OsiClpSolverInterface.hpp"
00026 #define OSI_SOLVER_INTERFACE OsiClpSolverInterface
00027
00028 using namespace std;
00029 using namespace flopc;
00030
00031
00033
00038 class StageNode {
00039 public:
00040 MP_set ASSETS;
00041 MP_index a;
00042 MP_variable x;
00043 MP_expression objFunction;
00044
00045 StageNode *ptParent;
00046 vector<StageNode *> children;
00047 int nodeNmb;
00048 double prob;
00049
00051
00056
00057 vector<VariableRef *> all_variables;
00058 vector<MP_constraint *> all_constraints;
00059
00060
00064 MP_constraint *pt_balance_constraint;
00066
00067 StageNode(StageNode *ptPred, const int nmbAssets, const int nodeN,
00068 const double condProb)
00069 : ASSETS(nmbAssets), x(ASSETS), ptParent(ptPred), nodeNmb(nodeN),
00070 prob(condProb)
00071 {
00072 if (ptParent != NULL) {
00073 prob *= ptParent->prob;
00074 ptParent->children.push_back(this);
00075 }
00076 }
00077 virtual ~StageNode() {
00078 for (int a = 0; a < (int) all_variables.size(); a++) {
00079 delete all_variables[a];
00080 }
00081 }
00082
00084
00090 virtual double get_wealth(const double *variableValues, const int nmbVars) {
00091 assert (nmbVars == (int) all_variables.size()
00092 && "Check that we have values of all variables");
00093 double wealth = 0;
00094 for (int i = 0; i < nmbVars; i++)
00095 wealth += variableValues[i];
00096 return wealth;
00097 }
00098
00099 protected:
00101
00103 virtual void make_obj_function_() {
00104 assert (children.size() > 0 && "Leaves should never call this function");
00105
00106 for (int i=0; i < (int) children.size(); i++) {
00107 children[i]->make_obj_function_();
00108 }
00109
00110
00111 objFunction = children[0]->objFunction;
00112 for (int i=1; i < (int) children.size(); i++) {
00113 objFunction = objFunction + children[i]->objFunction;
00114 }
00115 }
00116 };
00117
00119 class RootNode : public StageNode {
00120 public:
00121 MP_constraint initialBudget;
00122
00123 RootNode(const int nmbAssets, const double initWealth)
00124 : StageNode(NULL, nmbAssets, 0, 1.0)
00125 {
00126 initialBudget() = sum(ASSETS, x(ASSETS)) == initWealth;
00127
00128 for (int a = 0; a < nmbAssets; a++) {
00129 all_variables.push_back(new VariableRef(x(a)));
00130 }
00131 all_constraints.push_back(&initialBudget);
00132 pt_balance_constraint = NULL;
00133 }
00134
00136
00138 void make_objective_function() {
00139 make_obj_function_();
00140 }
00141 };
00142
00144 class MidStageNode : public StageNode {
00145 public:
00146 MP_constraint cashFlowBalance;
00147 MP_data Return;
00148
00150
00156 MidStageNode(StageNode *ptPred, const int nodeN, const double condProb,
00157 double *ptRetVect)
00158 : StageNode(ptPred, ptPred->ASSETS.size(), nodeN, condProb),
00159 Return(ptRetVect, ASSETS)
00160 {
00161
00162 cashFlowBalance = sum(ASSETS(a), ptParent->x(a) * Return(a))
00163 == sum(ASSETS(a), x(a));
00164
00165 for (int a = 0; a < ASSETS.size(); a++) {
00166 all_variables.push_back(new VariableRef(x(a)));
00167 }
00168 all_constraints.push_back(&cashFlowBalance);
00169 pt_balance_constraint = &cashFlowBalance;
00170 }
00171 };
00172
00174 class LeafNode : public StageNode {
00175 public:
00176 MP_variable w;
00177 MP_variable y;
00178 MP_constraint finalBalance;
00179 MP_data Return;
00180 double capTarget;
00181
00183
00187 LeafNode(StageNode *ptPred, const int nodeN, const double condProb,
00188 const double *ptRetVect, const double capTg)
00189 : StageNode(ptPred, ptPred->ASSETS.size(), nodeN, condProb), Return(ASSETS),
00190 capTarget(capTg)
00191 {
00192 Return.value(ptRetVect);
00193
00194 finalBalance = sum(ASSETS, ptParent->x(ASSETS) * Return(ASSETS))
00195 + w() - y() == capTarget;
00196
00197 all_variables.push_back(new VariableRef(w()));
00198 all_variables.push_back(new VariableRef(y()));
00199 all_constraints.push_back(&finalBalance);
00200 pt_balance_constraint = &finalBalance;
00201 }
00202
00204 double get_wealth(const double *variableValues, const int nmbVars) {
00205 assert (nmbVars == (int) all_variables.size()
00206 && "Check that we have values of all variables");
00207 return -variableValues[0] + variableValues[1] + capTarget;
00208 }
00209
00210 protected:
00212 void make_obj_function_() {
00213 objFunction = prob * (1.3 * w() - 1.1 * y());
00214 }
00215 };
00216
00217
00218
00219
00221 class ScenTreeStruct {
00222 public:
00223 int nmbNodes;
00224 int firstLeaf;
00225
00227 ScenTreeStruct(const int nNodes, const int firstL)
00228 : nmbNodes(nNodes), firstLeaf(firstL)
00229 {}
00230 virtual ~ScenTreeStruct(){}
00231
00233
00236 virtual int get_parent(int n) const = 0;
00237 };
00238
00240 class BinTreeStruct : public ScenTreeStruct {
00241 public:
00243 BinTreeStruct(const int nmbStages)
00244 : ScenTreeStruct(pow(2.0, nmbStages) - 1, pow(2.0, nmbStages-1) - 1) {}
00245
00246 int get_parent(int n) const {
00247 return (n-1) / 2;
00248 }
00249 };
00250
00251
00252
00253
00254 int main()
00255 {
00256
00257
00258
00259 const int nmbStages = 4;
00260 BinTreeStruct scTree(nmbStages);
00261
00262
00263 enum {stocks, bonds, nmbAssets};
00264 double InitBudget = 55;
00265 double CapTarget = 80;
00266
00267
00268 double retData[][nmbAssets] = {
00269 {1.25, 1.14},
00270 {1.06, 1.16},
00271 {1.21, 1.17},
00272 {1.07, 1.12},
00273 {1.15, 1.18},
00274 {1.06, 1.12},
00275 {1.26, 1.13},
00276 {1.07, 1.14},
00277 {1.25, 1.15},
00278 {1.06, 1.12},
00279 {1.05, 1.17},
00280 {1.06, 1.15},
00281 {1.05, 1.14},
00282 {1.06, 1.12}
00283 };
00284
00285
00286
00287
00288
00289
00290 MP_model &coreModel = MP_model::getDefaultModel();
00291 coreModel.setSolver(new OSI_SOLVER_INTERFACE);
00292 coreModel.silent();
00293
00294 int a, i, j, n, t;
00295
00296 vector<int> scenNodeNmb(nmbStages);
00297 n = scTree.firstLeaf;
00298 for (t = nmbStages; t > 0; t--) {
00299 scenNodeNmb[t-1] = n;
00300 n = scTree.get_parent(n);
00301 }
00302
00303
00304 vector<StageNode *> coreNodes(nmbStages);
00305 coreNodes[0] = new RootNode(nmbAssets, InitBudget);
00306 for (t = 1; t < nmbStages-1; t++) {
00307 coreNodes[t] = new MidStageNode(coreNodes[t-1], t, 1.0,
00308 retData[scenNodeNmb[t]-1]);
00309 }
00310 assert (t == nmbStages-1 && "t should be nmbStages-1 after the loop");
00311 coreNodes[t] = new LeafNode(coreNodes[t-1], t, 1.0, retData[scenNodeNmb[t]-1],
00312 CapTarget);
00313
00314
00315 RootNode *ptRoot = dynamic_cast<RootNode *>(coreNodes[0]);
00316
00317 ptRoot->make_objective_function();
00318 coreModel.setObjective(ptRoot->objFunction);
00319 coreModel.attach();
00320
00321
00322
00323 int nmbCoreCols = coreModel->getNumCols();
00324 int nmbCoreRows = coreModel->getNumRows();
00325
00326
00327 coreModel->writeMps("investment.core");
00328 cout << endl << "The core (deterministic) model has " << nmbCoreCols
00329 << " variables and " << nmbCoreRows << " constraints." << endl;
00330
00331
00332
00333 int *colStages = new int[nmbCoreCols];
00334 int checkSum = 0;
00335 for (t = 0; t < nmbStages; t++) {
00336 for (j = 0; j < (int) coreNodes[t]->all_variables.size(); j++) {
00337 int colIndx = coreNodes[t]->all_variables[j]->getColumn();
00338 #ifndef NDEBUG
00339 cout << "stage " << t << ": var no. " << j+1 << " is in column "
00340 << colIndx << endl;
00341 #endif
00342 colStages[colIndx] = t;
00343 checkSum += colIndx;
00344 }
00345 }
00346
00347 assert (checkSum == nmbCoreCols * (nmbCoreCols-1) / 2
00348 && "checkSum = sum of numbers from zero to nmbCoreCols-1");
00349
00350
00351 int *rowStages = new int[nmbCoreRows];
00352 checkSum = 0;
00353 for (t = 0; t < nmbStages; t++) {
00354 for (i = 0; i < (int) coreNodes[t]->all_constraints.size(); i++) {
00355 int rowIndx = *coreNodes[t]->all_constraints[i];
00356 #ifndef NDEBUG
00357 cout << "stage " << t << ": constraint no. " << i+1 << " is in row "
00358 << rowIndx << endl;
00359 #endif
00360 rowStages[rowIndx] = t;
00361 checkSum += rowIndx;
00362 }
00363 }
00364
00365 assert (checkSum == nmbCoreRows * (nmbCoreRows-1) / 2
00366 && "checkSum = sum of numbers from zero to nmbCoreRows-1");
00367
00368
00369 SmiCoreData stochCore(coreModel.operator->(), nmbStages,
00370 colStages, rowStages);
00371 delete[] colStages;
00372 delete[] rowStages;
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 SmiScnModel stochModel;
00384
00385
00386 CoinPackedMatrix ADiff;
00387
00388
00389 ADiff.reverseOrdering();
00390
00391 int nmbScen = scTree.nmbNodes - scTree.firstLeaf;
00392 double scenProb = 1.0 / nmbScen;
00393
00394
00395
00396
00397 for (int leaf = scTree.firstLeaf; leaf < scTree.nmbNodes; leaf++) {
00398 int scen = leaf - scTree.firstLeaf;
00399 n = leaf;
00400 t = nmbStages-1;
00401 cout << "Nodes in scenario " << scen + 1 << ": ";
00402 ADiff.clear();
00403 ADiff.setDimensions(nmbCoreRows, nmbCoreCols);
00404 while (n != scenNodeNmb[t]) {
00405 assert (n > 0 && t > 0 && "All scenarios must end in a common root");
00406 cout << setw(2) << n << " ";
00407 scenNodeNmb[t] = n;
00408
00409
00410 i = *coreNodes[t]->pt_balance_constraint;
00411 for (a=0; a<nmbAssets; a++) {
00412
00413 j = coreNodes[t-1]->x(a).getColumn();
00414 ADiff.modifyCoefficient(i, j, retData[n-1][a]);
00415 }
00416
00417
00418 n = scTree.get_parent(n);
00419 t--;
00420 }
00421 cout << endl;
00422
00423
00424 int branchStage = (scen == 0 ? 1 : t+1);
00425 int parentScen = (scen == 0 ? 0 : scen - 1);
00426 #ifdef NDEBUG
00427 stochModel.generateScenario(&stochCore, &ADiff,
00428 NULL, NULL, NULL, NULL, NULL,
00429 branchStage, parentScen, scenProb);
00430 #else
00431 int scenIndx = stochModel.generateScenario(&stochCore, &ADiff,
00432 NULL, NULL, NULL, NULL, NULL,
00433 branchStage, parentScen,
00434 scenProb);
00435 assert (scenIndx == scen && "Check index of the new scenario");
00436 #endif
00437 }
00438
00439
00440
00441
00442
00443
00444
00445 stochModel.setOsiSolverHandle(*new OSI_SOLVER_INTERFACE());
00446
00447
00448
00449 OsiSolverInterface *ptDetEqModel = stochModel.loadOsiSolverData();
00450
00451
00452 ptDetEqModel->writeMps("investment.det-equiv");
00453 cout << endl << "Solving the deterministic equivalent:" << endl;
00454 ptDetEqModel->initialSolve();
00455 cout << endl << "The deterministic equivalent model has "
00456 << ptDetEqModel->getNumCols() << " variables and "
00457 << ptDetEqModel->getNumRows() << " constraints." << endl;
00458 cout << "Optimal objective value = " << ptDetEqModel->getObjValue() << endl;
00459
00460
00461
00462
00463
00464 cout << endl << "The stochastic model has " << stochModel.getNumScenarios()
00465 << " scenarios." << endl;
00466 assert (stochModel.getNumScenarios() == nmbScen && "Check number of scens.");
00467
00468
00469 vector<double> nodeWealth(nmbStages, 0);
00470 double objValue = 0.0;
00471
00472
00473 const double *ptOsiSolution = ptDetEqModel->getColSolution();
00474
00475
00476 for (SmiScenarioIndex sc = 0; sc < nmbScen; sc ++) {
00477
00478 SmiScnNode *ptNode = stochModel.getLeafNode(sc);
00479 int nodeStage = nmbStages;
00480
00481 double scProb = ptNode->getModelProb();
00482 double scenObjVal = stochModel.getObjectiveValue(sc);
00483 objValue += scProb * scenObjVal;
00484 printf ("scen %d: prob = %.3f obj =%7.2f", sc+1, scProb, scenObjVal);
00485
00486
00487 while (ptNode != NULL) {
00488
00489 int startColInOsi = ptNode->getColStart();
00490 int nmbColsInOsi = ptNode->getNumCols();
00491
00492
00493 nodeWealth[nodeStage-1]
00494 = coreNodes[nodeStage-1]->get_wealth(&(ptOsiSolution[startColInOsi]),
00495 nmbColsInOsi);
00496
00497
00498 ptNode = ptNode->getParent();
00499 nodeStage--;
00500 }
00501
00502 printf ("; wealth:");
00503 for (int t = 0; t<nmbStages-1; t++)
00504 printf ("%6.2f ->", nodeWealth[t]);
00505 printf ("%6.2f\n", nodeWealth[nmbStages-1]);
00506 }
00507
00508 printf ("%15s Total obj = %7.3f\n", "", objValue);
00509 assert (fabs(objValue - ptDetEqModel->getObjValue()) < 1e-6 && "Check obj.");
00510
00511 for (t = 0; t < nmbStages; t++)
00512 delete coreNodes[t];
00513
00514 return 0;
00515 }