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
119
|
// Copyright (C) 2019 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#include "StOpt/regression/ContinuationCuts.h"
#include "StOpt/core/utils/comparisonUtils.h"
using namespace std;
using namespace Eigen;
namespace StOpt
{
ContinuationCuts::ContinuationCuts(const std::shared_ptr< SpaceGrid > &p_grid,
const std::shared_ptr< BaseRegression > &p_condExp,
const ArrayXXd &p_values) : m_grid(p_grid), m_condExp(p_condExp), m_regressedCutCoeff(p_grid->getDimension() + 1)
{
// nest on cuts
for (int ic = 0; ic < p_grid->getDimension() + 1; ++ic)
{
// coefficients of regressed functions (nb stock points, nb function basis)
// From Eigen 3.4 ArrayXXd valLoc = p_values(seqN(ic*p_condExp->getNbSimul(),p_condExp->getNbSimul()),all);
ArrayXXd valLoc = p_values.block(ic * p_condExp->getNbSimul(), 0, p_condExp->getNbSimul(), p_values.cols());
m_regressedCutCoeff(ic) = m_condExp->getCoordBasisFunctionMultiple(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> iterRegGrid = m_grid->getGridIterator();
while (iterRegGrid->isValid())
{
// coordinates
ArrayXd pointCoordReg = iterRegGrid->getCoordinate();
// point number
int ipoint = iterRegGrid->getCount();
// grid cuts
for (int id = 0 ; id < pointCoordReg.size(); ++id)
m_regressedCutCoeff(0).col(ipoint) -= m_regressedCutCoeff(id + 1).col(ipoint) * pointCoordReg(id);
iterRegGrid->next();
}
}
ArrayXXd ContinuationCuts::getCutsASim(const ArrayXXd &p_hypStock, const ArrayXd &p_coordinates) 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_condExp->getValue(p_coordinates, m_regressedCutCoeff(jc).col(ipoint));
}
iPointCut += 1;
}
iterRegGrid->next();
}
cuts.conservativeResize(nbCutsCoeff, iPointCut);
return cuts;
}
ArrayXXd ContinuationCuts::getCutsAllSimulations(const ArrayXXd &p_hypStock) const
{
int nbCutsCoeff = m_grid->getDimension() + 1;
// for return
ArrayXXd cuts(nbCutsCoeff * m_condExp->getNbSimul(), 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 * m_condExp->getNbSimul(), m_condExp->getNbSimul()) = m_condExp->reconstruction(m_regressedCutCoeff(jc).col(ipoint));
}
iPointCut += 1;
}
iterRegGrid->next();
}
cuts.conservativeResize(nbCutsCoeff * m_condExp->getNbSimul(), iPointCut);
return cuts;
}
}
|
