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
|
// @(#)root/mathmore:$Id$
// Author: L. Moneta June 2009
/**********************************************************************
* *
* Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
* *
* *
**********************************************************************/
// Implementation file for class MinimTransformFunction
#include "Math/MinimTransformFunction.h"
#include "Math/IFunctionfwd.h"
//#include <iostream>
#include <cmath>
#include <cassert>
namespace ROOT {
namespace Math {
MinimTransformFunction::MinimTransformFunction ( const IMultiGradFunction * f, const std::vector<EMinimVariableType> & types,
const std::vector<double> & values,
const std::map<unsigned int, std::pair<double, double> > & bounds) :
fX( values ),
fFunc(f)
{
// constructor of the class from a pointer to the function (which is managed)
// vector specifying the variable types (free, bounded or fixed, defined in enum EMinimVariableTypes )
// variable values (used for the fixed ones) and a map with the bounds (for the bounded variables)
unsigned int ntot = NTot(); // NTot is fFunc->NDim()
assert ( types.size() == ntot );
fVariables.reserve(ntot);
fIndex.reserve(ntot);
for (unsigned int i = 0; i < ntot; ++i ) {
if (types[i] == kFix )
fVariables.push_back( MinimTransformVariable( values[i]) );
else {
fIndex.push_back(i);
if ( types[i] == kDefault)
fVariables.push_back( MinimTransformVariable() );
else {
std::map<unsigned int, std::pair<double,double> >::const_iterator itr = bounds.find(i);
assert ( itr != bounds.end() );
double low = itr->second.first;
double up = itr->second.second;
if (types[i] == kBounds )
fVariables.push_back( MinimTransformVariable( low, up, new SinVariableTransformation()));
else if (types[i] == kLowBound)
fVariables.push_back( MinimTransformVariable( low, new SqrtLowVariableTransformation()));
else if (types[i] == kUpBound)
fVariables.push_back( MinimTransformVariable( up, new SqrtUpVariableTransformation()));
}
}
}
}
void MinimTransformFunction::Transformation( const double * x, double * xext) const {
// transform from internal to external
unsigned int nfree = fIndex.size();
// std::cout << "Transform: internal ";
// for (int i = 0; i < nfree; ++i) std::cout << x[i] << " ";
// std::cout << "\t\t";
for (unsigned int i = 0; i < nfree; ++i ) {
unsigned int extIndex = fIndex[i];
const MinimTransformVariable & var = fVariables[ extIndex ];
if (var.IsLimited() )
xext[ extIndex ] = var.InternalToExternal( x[i] );
else
xext[ extIndex ] = x[i];
}
// std::cout << "Transform: external ";
// for (int i = 0; i < fX.size(); ++i) std::cout << fX[i] << " ";
// std::cout << "\n";
}
void MinimTransformFunction::InvTransformation(const double * xExt, double * xInt) const {
// inverse function transformation (external -> internal)
for (unsigned int i = 0; i < NDim(); ++i ) {
unsigned int extIndex = fIndex[i];
const MinimTransformVariable & var = fVariables[ extIndex ];
assert ( !var.IsFixed() );
if (var.IsLimited() )
xInt[ i ] = var.ExternalToInternal( xExt[extIndex] );
else
xInt[ i ] = xExt[extIndex];
}
}
void MinimTransformFunction::InvStepTransformation(const double * x, const double * sExt, double * sInt) const {
// inverse function transformation for steps (external -> internal)
for (unsigned int i = 0; i < NDim(); ++i ) {
unsigned int extIndex = fIndex[i];
const MinimTransformVariable & var = fVariables[ extIndex ];
assert ( !var.IsFixed() );
if (var.IsLimited() ) {
// bound variables
double x2 = x[extIndex] + sExt[extIndex];
if (var.HasUpperBound() && x2 >= var.UpperBound() )
x2 = x[extIndex] - sExt[extIndex];
// transform x and x2
double xint = var.ExternalToInternal ( x[extIndex] );
double x2int = var.ExternalToInternal( x2 );
sInt[i] = std::abs( x2int - xint);
}
else
sInt[ i ] = sExt[extIndex];
}
}
void MinimTransformFunction::GradientTransformation(const double * x, const double *gExt, double * gInt) const {
//transform gradient vector (external -> internal) at internal point x
unsigned int nfree = fIndex.size();
for (unsigned int i = 0; i < nfree; ++i ) {
unsigned int extIndex = fIndex[i];
const MinimTransformVariable & var = fVariables[ extIndex ];
assert (!var.IsFixed() );
if (var.IsLimited() )
gInt[i] = gExt[ extIndex ] * var.DerivativeIntToExt( x[i] );
else
gInt[i] = gExt[ extIndex ];
}
}
void MinimTransformFunction::MatrixTransformation(const double * x, const double *covInt, double * covExt) const {
//transform covariance matrix (internal -> external) at internal point x
// use row storages for matrices m(i,j) = rep[ i * dim + j]
// ignore fixed points
unsigned int nfree = fIndex.size();
unsigned int ntot = NTot();
for (unsigned int i = 0; i < nfree; ++i ) {
unsigned int iext = fIndex[i];
const MinimTransformVariable & ivar = fVariables[ iext ];
assert (!ivar.IsFixed());
double ddi = ( ivar.IsLimited() ) ? ivar.DerivativeIntToExt( x[i] ) : 1.0;
// loop on j variables for not fixed i variables (forget that matrix is symmetric) - could be optimized
for (unsigned int j = 0; j < nfree; ++j ) {
unsigned int jext = fIndex[j];
const MinimTransformVariable & jvar = fVariables[ jext ];
double ddj = ( jvar.IsLimited() ) ? jvar.DerivativeIntToExt( x[j] ) : 1.0;
assert (!jvar.IsFixed() );
covExt[ iext * ntot + jext] = ddi * ddj * covInt[ i * nfree + j];
}
}
}
} // end namespace Math
} // end namespace ROOT
|
