2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4 // + This file is part of enGrid. +
6 // + Copyright 2008-2010 enGits GmbH +
8 // + enGrid is free software: you can redistribute it and/or modify +
9 // + it under the terms of the GNU General Public License as published by +
10 // + the Free Software Foundation, either version 3 of the License, or +
11 // + (at your option) any later version. +
13 // + enGrid is distributed in the hope that it will be useful, +
14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
16 // + GNU General Public License for more details. +
18 // + You should have received a copy of the GNU General Public License +
19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23 #include "optimisation.h"
24 #include "guimainwindow.h"
26 ErrorFunction::ErrorFunction()
37 void ErrorFunction::set(QString settings_txt
)
39 QStringList items
= settings_txt
.split(';');
40 if (items
.size() != 2) {
41 EG_ERR_RETURN("syntax error for error weighting");
43 m_Err0
= items
[0].trimmed().toDouble();
45 m_XS
= items
[1].trimmed().toDouble();
48 m_Exp
= log(m_ErrS
/m_Err0
)/log(1.0 - m_XS
);
52 double ErrorFunction::operator()(double x
)
56 double x0
= fabs(1-x
);
57 e
= m_Err0
*pow(x0
, m_Exp
);
58 m_MaxErr
= max(m_MaxErr
, x0
);
67 double ErrorFunction::averageError()
69 if (m_NumTotalCalls
== 0) {
72 return m_AverageError
/m_NumTotalCalls
;
75 double ErrorFunction::totalError()
77 if (m_NumCalls
== 0) {
80 return m_TotalError
/m_NumCalls
;
83 void ErrorFunction::reset(bool reset_average
)
94 Optimisation::Optimisation()
98 for (int i
= 0; i
< 3; ++i
) {
99 F
[i
] = new double*[3];
100 for (int j
= 0; j
< 3; ++j
) {
101 F
[i
][j
] = new double[3];
104 m_ErrorFunctions
.clear();
107 void Optimisation::getErrSet(QString group
, QString key
, double err0
, double xs
, ErrorFunction
* err_func
)
109 QString err0_txt
, xs_txt
;
110 err0_txt
.setNum(err0
);
112 QString value
= err0_txt
+ "; " + xs_txt
;
114 QSettings
*qset
= GuiMainWindow::settings();
115 QString typed_key
= "string/" + key
;
116 if (group
!= QObject::tr("General")) {
117 qset
->beginGroup(group
);
119 if (!qset
->contains(typed_key
)) {
120 qset
->setValue(typed_key
, value
);
122 variable
= (qset
->value(typed_key
,variable
)).toString();
123 if (group
!= QObject::tr("General")) {
126 err_func
->set(variable
);
127 err_func
->setName(key
.replace(" ", "_"));
128 m_ErrorFunctions
.append(err_func
);
131 double Optimisation::angleX(const vec3_t
&v1
, const vec3_t
&v2
)
133 double va1
= v1
.abs();
134 double va2
= v2
.abs();
135 double scal
= max(0.0, v1
*v2
);
136 double alpha
= acos(scal
/(va1
*va2
));
137 return fabs(1 - 2*alpha
/M_PI
);
140 void Optimisation::computeDerivatives(vec3_t x
)
142 F
[0][0][0] = func(x
[0]-Dx
, x
[1]-Dy
, x
[2]-Dz
);
143 F
[1][0][0] = func(x
[0], x
[1]-Dy
, x
[2]-Dz
);
144 F
[2][0][0] = func(x
[0]+Dx
, x
[1]-Dy
, x
[2]-Dz
);
145 F
[0][1][0] = func(x
[0]-Dx
, x
[1], x
[2]-Dz
);
146 F
[1][1][0] = func(x
[0], x
[1], x
[2]-Dz
);
147 F
[2][1][0] = func(x
[0]+Dx
, x
[1], x
[2]-Dz
);
148 F
[0][2][0] = func(x
[0]-Dx
, x
[1]+Dy
, x
[2]-Dz
);
149 F
[1][2][0] = func(x
[0], x
[1]+Dy
, x
[2]-Dz
);
150 F
[2][2][0] = func(x
[0]+Dx
, x
[1]+Dy
, x
[2]-Dz
);
151 F
[0][0][1] = func(x
[0]-Dx
, x
[1]-Dy
, x
[2]);
152 F
[1][0][1] = func(x
[0], x
[1]-Dy
, x
[2]);
153 F
[2][0][1] = func(x
[0]+Dx
, x
[1]-Dy
, x
[2]);
154 F
[0][1][1] = func(x
[0]-Dx
, x
[1], x
[2]);
155 F
[1][1][1] = func(x
[0], x
[1], x
[2]);
156 F
[2][1][1] = func(x
[0]+Dx
, x
[1], x
[2]);
157 F
[0][2][1] = func(x
[0]-Dx
, x
[1]+Dy
, x
[2]);
158 F
[1][2][1] = func(x
[0], x
[1]+Dy
, x
[2]);
159 F
[2][2][1] = func(x
[0]+Dx
, x
[1]+Dy
, x
[2]);
160 F
[0][0][2] = func(x
[0]-Dx
, x
[1]-Dy
, x
[2]+Dz
);
161 F
[1][0][2] = func(x
[0], x
[1]-Dy
, x
[2]+Dz
);
162 F
[2][0][2] = func(x
[0]+Dx
, x
[1]-Dy
, x
[2]+Dz
);
163 F
[0][1][2] = func(x
[0]-Dx
, x
[1], x
[2]+Dz
);
164 F
[1][1][2] = func(x
[0], x
[1], x
[2]+Dz
);
165 F
[2][1][2] = func(x
[0]+Dx
, x
[1], x
[2]+Dz
);
166 F
[0][2][2] = func(x
[0]-Dx
, x
[1]+Dy
, x
[2]+Dz
);
167 F
[1][2][2] = func(x
[0], x
[1]+Dy
, x
[2]+Dz
);
168 F
[2][2][2] = func(x
[0]+Dx
, x
[1]+Dy
, x
[2]+Dz
);
170 grad_f
[0] = (F
[2][1][1]-F
[0][1][1])/(2*Dx
);
171 grad_f
[1] = (F
[1][2][1]-F
[1][0][1])/(2*Dy
);
172 grad_f
[2] = (F
[1][1][2]-F
[1][1][0])/(2*Dz
);
174 J
[0][0] = (F
[0][1][1]-2*F
[1][1][1]+F
[2][1][1])/(Dx
*Dx
);
175 J
[1][1] = (F
[1][0][1]-2*F
[1][1][1]+F
[1][2][1])/(Dy
*Dy
);
176 J
[2][2] = (F
[1][1][0]-2*F
[1][1][1]+F
[1][1][2])/(Dz
*Dz
);
178 J
[0][1] = ((F
[2][2][1]-F
[0][2][1]) - (F
[2][0][1]-F
[0][0][1]))/(4*Dx
*Dy
);
179 J
[0][2] = ((F
[2][1][2]-F
[0][1][2]) - (F
[2][1][0]-F
[0][1][0]))/(4*Dx
*Dz
);
180 J
[1][2] = ((F
[1][2][2]-F
[1][2][2]) - (F
[1][0][0]-F
[1][0][0]))/(4*Dy
*Dz
);
186 vec3_t
Optimisation::optimise(vec3_t x
)
188 computeDerivatives(x
);
189 mat3_t M
= J
.inverse();
190 return (-1)*(M
*grad_f
);
193 void Optimisation::resetErrorFunctions(bool reset_average
)
195 foreach (ErrorFunction
*err_func
, m_ErrorFunctions
) {
196 err_func
->reset(reset_average
);
200 double Optimisation::totalError()
203 foreach (ErrorFunction
*err_func
, m_ErrorFunctions
) {
204 e
+= err_func
->totalError();
209 void Optimisation::printErrors()
211 foreach (ErrorFunction
*err_func
, m_ErrorFunctions
) {
212 if (err_func
->active()) {
213 cout
<< qPrintable(err_func
->name()) << ":\n average=" << err_func
->averageError() << "\n maximum=" << err_func
->maxError() << endl
;