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
|
// Copyright (C) 2017 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#include "StOpt/regression/BaseRegression.h"
using namespace std;
using namespace Eigen;
namespace StOpt
{
BaseRegression::BaseRegression(): m_bRotationAndRescale(false) {}
BaseRegression::BaseRegression(const bool &p_bRotationAndRescale):
m_bRotationAndRescale(p_bRotationAndRescale) {}
BaseRegression::BaseRegression(const bool &p_bZeroDate,
const ArrayXXd &p_particles,
const bool &p_bRotationAndRescale) :
m_bZeroDate(p_bZeroDate), m_bRotationAndRescale(p_bRotationAndRescale),
m_meanX(p_particles.rows()),
m_etypX(p_particles.rows()),
m_particles(p_particles)
{
preProcessData();
}
BaseRegression::BaseRegression(const bool &p_bZeroDate, const bool &p_bRotationAndRescale) :
m_bZeroDate(p_bZeroDate), m_bRotationAndRescale(p_bRotationAndRescale)
{}
BaseRegression:: BaseRegression(const bool &p_bZeroDate,
const ArrayXd &p_meanX,
const ArrayXd &p_etypX,
const MatrixXd &p_svdMatrix,
const bool &p_bRotationAndRescale) :
m_bZeroDate(p_bZeroDate), m_bRotationAndRescale(p_bRotationAndRescale)
{
if (p_meanX.size() > 0)
{
m_meanX = p_meanX;
m_etypX = p_etypX;
m_svdMatrix = p_svdMatrix;
}
}
BaseRegression::BaseRegression(const BaseRegression &p_object):
m_bZeroDate(p_object.getBZeroDate()),
m_bRotationAndRescale(p_object.getBRotationAndRescale()),
m_meanX(p_object.getMeanX()),
m_etypX(p_object.getEtypX()),
m_svdMatrix(p_object.getSvdMatrix()),
m_sing(p_object.getSing()),
m_particles(p_object.getParticles())
{}
void BaseRegression::preProcessData()
{
m_meanX = ArrayXd::Zero(m_particles.rows());
m_etypX = ArrayXd::Constant(m_particles.rows(), 1.);
m_svdMatrix = MatrixXd::Identity(m_particles.rows(), m_particles.rows());
if ((m_bRotationAndRescale) && (!m_bZeroDate))
{
// renormalize data
for (int id = 0; id < m_particles.rows(); ++id)
{
m_meanX(id) = m_particles.row(id).mean();
m_particles.row(id) -= m_meanX(id);
m_etypX(id) = sqrt(m_particles.row(id).pow(2).mean());
m_particles.row(id) /= m_etypX(id);
}
BDCSVD<MatrixXd> svd(m_particles.matrix(), ComputeThinU);
m_sing = svd.singularValues().array();
MatrixXd svdMatrix = svd.matrixU();
for (int id = 0; id < m_particles.rows(); ++id)
{
double vNrom = svdMatrix(id, id) / fabs(svdMatrix(id, id));
for (int iid = 0; iid < m_particles.rows(); ++iid)
svdMatrix(iid, id) *= vNrom;
}
// transpose
m_svdMatrix = svdMatrix.transpose();
m_particles = (m_svdMatrix * m_particles.matrix()).array(); // rotation
}
}
void BaseRegression::updateSimulationsBase(const bool &p_bZeroDate,
const ArrayXXd &p_particles)
{
m_bZeroDate = p_bZeroDate;
m_particles = p_particles;
preProcessData();
}
ArrayXd BaseRegression:: getParticle(const int &p_iPart) const
{
assert((p_iPart == 0) || (p_iPart < m_particles.cols()));
if ((p_iPart == 0) && (m_particles.size() == 0))
return ArrayXd();
else
return m_particles.col(p_iPart);
}
}
|
