1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
// Copyright (C) 2019 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#include "StOpt/tree/ContinuationCutsTree.h"
#include "StOpt/core/utils/comparisonUtils.h"
using namespace std;
using namespace Eigen;
namespace StOpt
{
ContinuationCutsTree::ContinuationCutsTree(const std::shared_ptr< SpaceGrid > &p_grid,
const std::shared_ptr< Tree > &p_condExp,
const ArrayXXd &p_values) : m_grid(p_grid), m_cutCoeff(p_grid->getDimension() + 1)
{
// nest on cuts
for (int ic = 0; ic < p_grid->getDimension() + 1; ++ic)
{
ArrayXXd valLoc = p_values.block(ic * p_condExp->getNbNodesNextDate(), 0, p_condExp->getNbNodesNextDate(), p_values.cols());
m_cutCoeff[ic] = p_condExp->expCondMultiple(valLoc.transpose()).transpose();
}
// for first coefficient cuts calculate \f$ \bar a_0 = a_0 - \sum_{i=1}^d a_i \bar x_i \f$
// iterator
shared_ptr<GridIterator> iterGrid = m_grid->getGridIterator();
while (iterGrid->isValid())
{
// coordinates
ArrayXd pointCoord = iterGrid->getCoordinate();
// point number
int ipoint = iterGrid->getCount();
// grid cuts
for (int id = 0 ; id < pointCoord.size(); ++id)
m_cutCoeff[0].col(ipoint) -= m_cutCoeff[id + 1].col(ipoint) * pointCoord(id);
iterGrid->next();
}
}
ArrayXXd ContinuationCutsTree::getCutsANode(const ArrayXXd &p_hypStock, const int &p_node) const
{
int nbCutsCoeff = m_grid->getDimension() + 1;
// for return
ArrayXXd cuts(nbCutsCoeff, m_grid->getNbPoints());
int iPointCut = 0;
// nest on grid points
shared_ptr<GridIterator> iterRegGrid = m_grid->getGridIterator();
while (iterRegGrid->isValid())
{
// coordinates
Eigen::ArrayXd pointCoordReg = iterRegGrid->getCoordinate();
// test if inside the hypercube
bool bInside = true;
for (int id = 0 ; id < pointCoordReg.size(); ++id)
if (isStrictlyLesser(pointCoordReg(id), p_hypStock(id, 0)) || (isStrictlyMore(pointCoordReg(id), p_hypStock(id, 1))))
{
bInside = false;
break;
}
if (bInside)
{
// point number
int ipoint = iterRegGrid->getCount();
for (int jc = 0; jc < nbCutsCoeff; ++jc)
{
// reconstruct the value for all simulations
cuts(jc, iPointCut) = m_cutCoeff[jc](p_node, ipoint);
}
iPointCut += 1;
}
iterRegGrid->next();
}
cuts.conservativeResize(nbCutsCoeff, iPointCut);
return cuts;
}
ArrayXXd ContinuationCutsTree::getCutsAllNodes(const ArrayXXd &p_hypStock) const
{
int nbCutsCoeff = m_grid->getDimension() + 1;
int nbNodes = m_cutCoeff[0].rows();
// for return
ArrayXXd cuts(nbCutsCoeff * nbNodes, m_grid->getNbPoints());
int iPointCut = 0;
// nest on grid points
shared_ptr<GridIterator> iterRegGrid = m_grid->getGridIterator();
while (iterRegGrid->isValid())
{
// coordinates
Eigen::ArrayXd pointCoordReg = iterRegGrid->getCoordinate();
// test if inside the hypercube
bool bInside = true;
for (int id = 0 ; id < pointCoordReg.size(); ++id)
if (isStrictlyLesser(pointCoordReg(id), p_hypStock(id, 0)) || (isStrictlyMore(pointCoordReg(id), p_hypStock(id, 1))))
{
bInside = false;
break;
}
if (bInside)
{
// point number
int ipoint = iterRegGrid->getCount();
for (int jc = 0; jc < nbCutsCoeff; ++jc)
{
// reconstruct the value for all simulations
cuts.col(iPointCut).segment(jc * nbNodes, nbNodes) = m_cutCoeff[jc].col(ipoint);
}
iPointCut += 1;
}
iterRegGrid->next();
}
cuts.conservativeResize(nbCutsCoeff * nbNodes, iPointCut);
return cuts;
}
}
|
