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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
|
// Copyright (C) 2020 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#include <vector>
#include <memory>
#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen ;
namespace StOpt
{
/// \brief Calculate local sums as explained in Langrené-Warin 2020
/// \param p_x particules (sample) size : (dimension, nbSim)
/// \param p_z rectilinear points
/// \param p_y estimate for each p_x point
/// \param p_nbPtZ helper to store index in local sum
/// \param p_nbPtZT total number of local sum
ArrayXd localSumCalc(const ArrayXXd &p_x,
const vector< shared_ptr<ArrayXd> > &p_z,
const ArrayXd &p_y,
const ArrayXi &p_nbPtZ, const int &p_nbPtZT)
{
ArrayXXi iSort(p_x.rows(), p_x.cols());
ArrayXXd sX(p_x.rows(), p_x.cols());
for (int id = 0 ; id < p_x.rows(); ++id)
{
vector<pair<double, int> > toSort(p_x.cols());
for (int is = 0; is < p_x.cols(); ++is)
toSort[is] = make_pair(p_x(id, is), is);
sort(toSort.begin(), toSort.end());
for (int is = 0; is < p_x.cols(); ++is)
{
iSort(id, is) = toSort[is].second;
sX(id, is) = toSort[is].first;
}
}
ArrayXXi idx(p_x.rows(), p_x.cols());
// loop on dimension
for (int id = 0; id < p_x.rows(); ++id)
{
int xidx = 0;
int zidx = 0;
// select only important
while ((xidx < p_x.cols()) && (sX(id, xidx) <= (*p_z[id])(p_z[id]->size() - 1)))
{
if (sX(id, xidx) <= (*p_z[id])(zidx))
{
idx(id, iSort(id, xidx)) = zidx;
xidx += 1;
}
else
{
zidx += 1;
}
}
for (int ix = xidx; ix < p_x.cols(); ++ix)
{
idx(id, iSort(id, ix)) = -1;
}
}
// return
ArrayXd localSum = ArrayXd::Zero(p_nbPtZT);
for (int is = 0; is < p_x.cols(); ++is)
{
// index
if (idx.col(is).minCoeff() >= 0)
{
int iindex = idx(0, is) ;
for (int id = 1; id < p_x.rows() ; ++id)
iindex += idx(id, is) * p_nbPtZ(id);
localSum(iindex) += p_y(is);
}
}
localSum /= p_x.cols();
return localSum;
}
/// \brief Calculate fast CDF using Langrené Warin 2020 method
/// \param p_x particules (sample) size : (dimension, nbSim)
/// \param p_z rectilinear points
/// \param p_y estimate for each p_x point
/// \return for each point of the grid return the CDF
ArrayXd fastCDF(const ArrayXXd &p_x, const vector< shared_ptr<ArrayXd> > &p_z, const ArrayXd &p_y)
{
// store nbpt per dimension before
ArrayXi nbPtZ(p_z.size() + 1);
nbPtZ(0) = 1;
for (size_t id = 0; id < p_z.size(); ++id)
nbPtZ(id + 1) = nbPtZ(id) * p_z[id]->size();
// nb point
int nbPtZT = nbPtZ(p_z.size());
// helper
ArrayXi nbPtZInv(p_z.size());
nbPtZInv(p_z.size() - 1) = 1;
for (int id = p_z.size() - 2; id >= 0 ; --id)
nbPtZInv(id) = nbPtZInv(id + 1) * p_z[id + 1]->size();
// calculate local sum
ArrayXd localSum = localSumCalc(p_x, p_z, p_y, nbPtZ, nbPtZT);
// store index
ArrayXi index = ArrayXi::Constant(p_z.size(), -1);
// for return
ArrayXd retCDF = ArrayXd::Zero(nbPtZT);
// to store partial sum
vector< shared_ptr< ArrayXd> > sumOfLocal(p_z.size() + 1);
int isumOfLocal = nbPtZT;
sumOfLocal[0] = make_shared< ArrayXd>(ArrayXd::Zero(isumOfLocal));
for (size_t id = 1; id <= p_z.size(); ++id)
{
isumOfLocal /= p_z[id - 1]->size();
sumOfLocal[id] = make_shared< ArrayXd>(ArrayXd::Zero(isumOfLocal));
}
(*sumOfLocal[0]) = localSum;
// previous index
ArrayXi indexPrev(p_z.size());
// nest on number of Z points
for (int iloc = 0; iloc < nbPtZT; ++iloc)
{
indexPrev = index;
// index calculation for CDF
int ipt = iloc;
for (size_t id = 0; id < p_z.size(); ++id)
{
index(id) = static_cast<int>(ipt / nbPtZInv(id));
ipt -= index(id) * nbPtZInv(id);
}
for (size_t id = 0; id < p_z.size(); ++id)
{
if (index(id) != indexPrev(id))
{
// update sum of local
for (int ipt = 0; ipt < sumOfLocal[id + 1]->size(); ++ipt)
{
int iindexInSum = index(id) + ipt * p_z[id]->size();
// local index in sumOfLocal array
(*sumOfLocal[id + 1])(ipt) += (*sumOfLocal[id])(iindexInSum);
}
// update sum of local in other dimensions
for (size_t idN = id + 1; idN < p_z.size(); ++idN)
{
// loop on points
for (int ipt = 0; ipt < sumOfLocal[idN + 1]->size(); ++ipt)
{
(*sumOfLocal[idN + 1])(ipt) = (*sumOfLocal[idN])(ipt * p_z[idN]->size());
}
}
break;
}
}
// position in return array
int illoc = index(0);
for (size_t id = 1; id < p_z.size(); ++id)
illoc += index(id) * nbPtZ(id);
// store result
retCDF(illoc) = (*sumOfLocal[p_z.size()])(0);
}
return retCDF;
}
}
|
