src/libengrid/optimisation.cpp

   1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   2 // +                                                                      +
   3 // + This file is part of enGrid.                                         +
   4 // +                                                                      +
   5 // + Copyright 2008-2014 enGits GmbH                                      +
   6 // +                                                                      +
   7 // + enGrid is free software: you can redistribute it and/or modify       +
   8 // + it under the terms of the GNU General Public License as published by +
   9 // + the Free Software Foundation, either version 3 of the License, or    +
  10 // + (at your option) any later version.                                  +
  11 // +                                                                      +
  12 // + enGrid is distributed in the hope that it will be useful,            +
  13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  15 // + GNU General Public License for more details.                         +
  16 // +                                                                      +
  17 // + You should have received a copy of the GNU General Public License    +
  18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  19 // +                                                                      +
  20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  21 #include "optimisation.h"
  22 #include "guimainwindow.h"
  23
  24 ErrorFunction::ErrorFunction()
  25 {
  26   m_Err0 = 1.0;
  27   m_ErrS = 0.5;
  28   m_XS   = 0.5;
  29   m_Exp  = 1.0;
  30
  31   m_TotalError = 0.0;
  32   m_Active     = true;
  33 }
  34
  35 void ErrorFunction::set(QString settings_txt)
  36 {
  37   QStringList items = settings_txt.split(';');
  38   if (items.size() != 2) {
  39     EG_ERR_RETURN("syntax error for error weighting");
  40   }
  41   m_Err0 = items[0].trimmed().toDouble();
  42   m_ErrS = 0.5*m_Err0;
  43   m_XS   = items[1].trimmed().toDouble();
  44   m_Exp = 1.0;
  45   if (m_Err0 > 1e-3) {
  46     m_Exp = log(m_ErrS/m_Err0)/log(1.0 - m_XS);
  47   }
  48 }
  49
  50 double ErrorFunction::operator()(double x)
  51 {
  52   double e = 0;
  53   if (m_Active) {
  54     double x0 = fabs(1-x);
  55     e = m_Err0*pow(x0, m_Exp);
  56     m_MaxErr = max(m_MaxErr, x0);
  57     m_TotalError += e;
  58     m_AverageError += x0;
  59     ++m_NumCalls;
  60     ++m_NumTotalCalls;
  61   }
  62   return e;
  63 }
  64
  65 double ErrorFunction::averageError()
  66 {
  67   if (m_NumTotalCalls == 0) {
  68     return 0;
  69   }
  70   return m_AverageError/m_NumTotalCalls;
  71 }
  72
  73 double ErrorFunction::totalError()
  74 {
  75   if (m_NumCalls == 0) {
  76     return 0;
  77   }
  78   return m_TotalError/m_NumCalls;
  79 }
  80
  81 void ErrorFunction::reset(bool reset_average)
  82 {
  83   m_MaxErr = 0;
  84   m_TotalError = 0;
  85   m_NumCalls = 0;
  86   if (reset_average) {
  87     m_AverageError = 0;
  88     m_NumTotalCalls = 0;
  89   }
  90 }
  91
  92 Optimisation::Optimisation()
  93 {
  94   setDeltas(1e-6);
  95   F = new double**[3];
  96   for (int i = 0; i < 3; ++i) {
  97     F[i] = new double*[3];
  98     for (int j = 0; j < 3; ++j) {
  99       F[i][j] = new double[3];
 100     }
 101   }
 102   m_ErrorFunctions.clear();
 103 }
 104
 105 void Optimisation::getErrSet(QString group, QString key, double err0, double xs, ErrorFunction* err_func)
 106 {
 107   QString err0_txt, xs_txt;
 108   err0_txt.setNum(err0);
 109   xs_txt.setNum(xs);
 110   QString value = err0_txt + "; " + xs_txt;
 111   QString variable;
 112   QSettings *qset = GuiMainWindow::settings();
 113   QString typed_key = "string/" + key;
 114   if (group != QObject::tr("General")) {
 115     qset->beginGroup(group);
 116   }
 117   if (!qset->contains(typed_key)) {
 118     qset->setValue(typed_key, value);
 119   }
 120   variable = (qset->value(typed_key,variable)).toString();
 121   if (group != QObject::tr("General")) {
 122     qset->endGroup();
 123   }
 124   err_func->set(variable);
 125   err_func->setName(key.replace(" ", "_"));
 126   m_ErrorFunctions.append(err_func);
 127 }
 128
 129 double Optimisation::angleX(const vec3_t &v1, const vec3_t &v2)
 130 {
 131   double va1  = v1.abs();
 132   double va2  = v2.abs();
 133   double scal = max(0.0, v1*v2);
 134   double alpha = acos(scal/(va1*va2));
 135   return fabs(1 - 2*alpha/M_PI);
 136 }
 137
 138 void Optimisation::computeDerivatives(vec3_t x)
 139 {
 140   F[0][0][0] = func(x[0]-Dx, x[1]-Dy, x[2]-Dz);
 141   F[1][0][0] = func(x[0],    x[1]-Dy, x[2]-Dz);
 142   F[2][0][0] = func(x[0]+Dx, x[1]-Dy, x[2]-Dz);
 143   F[0][1][0] = func(x[0]-Dx, x[1],    x[2]-Dz);
 144   F[1][1][0] = func(x[0],    x[1],    x[2]-Dz);
 145   F[2][1][0] = func(x[0]+Dx, x[1],    x[2]-Dz);
 146   F[0][2][0] = func(x[0]-Dx, x[1]+Dy, x[2]-Dz);
 147   F[1][2][0] = func(x[0],    x[1]+Dy, x[2]-Dz);
 148   F[2][2][0] = func(x[0]+Dx, x[1]+Dy, x[2]-Dz);
 149   F[0][0][1] = func(x[0]-Dx, x[1]-Dy, x[2]);
 150   F[1][0][1] = func(x[0],    x[1]-Dy, x[2]);
 151   F[2][0][1] = func(x[0]+Dx, x[1]-Dy, x[2]);
 152   F[0][1][1] = func(x[0]-Dx, x[1],    x[2]);
 153   F[1][1][1] = func(x[0],    x[1],    x[2]);
 154   F[2][1][1] = func(x[0]+Dx, x[1],    x[2]);
 155   F[0][2][1] = func(x[0]-Dx, x[1]+Dy, x[2]);
 156   F[1][2][1] = func(x[0],    x[1]+Dy, x[2]);
 157   F[2][2][1] = func(x[0]+Dx, x[1]+Dy, x[2]);
 158   F[0][0][2] = func(x[0]-Dx, x[1]-Dy, x[2]+Dz);
 159   F[1][0][2] = func(x[0],    x[1]-Dy, x[2]+Dz);
 160   F[2][0][2] = func(x[0]+Dx, x[1]-Dy, x[2]+Dz);
 161   F[0][1][2] = func(x[0]-Dx, x[1],    x[2]+Dz);
 162   F[1][1][2] = func(x[0],    x[1],    x[2]+Dz);
 163   F[2][1][2] = func(x[0]+Dx, x[1],    x[2]+Dz);
 164   F[0][2][2] = func(x[0]-Dx, x[1]+Dy, x[2]+Dz);
 165   F[1][2][2] = func(x[0],    x[1]+Dy, x[2]+Dz);
 166   F[2][2][2] = func(x[0]+Dx, x[1]+Dy, x[2]+Dz);
 167
 168   grad_f[0] = (F[2][1][1]-F[0][1][1])/(2*Dx);
 169   grad_f[1] = (F[1][2][1]-F[1][0][1])/(2*Dy);
 170   grad_f[2] = (F[1][1][2]-F[1][1][0])/(2*Dz);
 171
 172   J[0][0] = (F[0][1][1]-2*F[1][1][1]+F[2][1][1])/(Dx*Dx);
 173   J[1][1] = (F[1][0][1]-2*F[1][1][1]+F[1][2][1])/(Dy*Dy);
 174   J[2][2] = (F[1][1][0]-2*F[1][1][1]+F[1][1][2])/(Dz*Dz);
 175
 176   J[0][1] = ((F[2][2][1]-F[0][2][1]) - (F[2][0][1]-F[0][0][1]))/(4*Dx*Dy);
 177   J[0][2] = ((F[2][1][2]-F[0][1][2]) - (F[2][1][0]-F[0][1][0]))/(4*Dx*Dz);
 178   J[1][2] = ((F[1][2][2]-F[1][2][2]) - (F[1][0][0]-F[1][0][0]))/(4*Dy*Dz);
 179   J[1][0] = J[0][1];
 180   J[2][0] = J[0][2];
 181   J[2][1] = J[1][2];
 182 }
 183
 184 vec3_t Optimisation::optimise(vec3_t x)
 185 {
 186   computeDerivatives(x);
 187   mat3_t M = J.inverse();
 188   vec3_t x_new = (-1)*(M*grad_f);
 189   if (!checkVector(x_new)) {
 190     x_new = x;
 191   }
 192   return x;
 193 }
 194
 195
 196 void Optimisation::resetErrorFunctions(bool reset_average)
 197 {
 198   foreach (ErrorFunction *err_func, m_ErrorFunctions) {
 199     err_func->reset(reset_average);
 200   }
 201 }
 202
 203 double Optimisation::totalError()
 204 {
 205   double e = 0;
 206   foreach (ErrorFunction *err_func, m_ErrorFunctions) {
 207     e += err_func->totalError();
 208   }
 209   return e;
 210 }
 211
 212 void Optimisation::printErrors()
 213 {
 214   foreach (ErrorFunction *err_func, m_ErrorFunctions) {
 215     if (err_func->active()) {
 216       cout << qPrintable(err_func->name()) << ":\n  average=" << err_func->averageError() << "\n  maximum=" << err_func->maxError() << endl;
 217     }
 218   }
 219 }