src/optimisation.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008,2009 Oliver Gloth                                     +
   7 // +                                                                      +
   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.                                  +
  12 // +                                                                      +
  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.                         +
  17 // +                                                                      +
  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/>.       +
  20 // +                                                                      +
  21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 //
  23 #include "optimisation.h"
  24 #include "guimainwindow.h"
  25
  26 ErrorFunction::ErrorFunction()
  27 {
  28   m_Err0 = 1.0;
  29   m_ErrS = 0.5;
  30   m_XS   = 0.5;
  31   m_Exp  = 1.0;
  32
  33   m_TotalError = 0.0;
  34   m_Active     = true;
  35 }
  36
  37 void ErrorFunction::set(QString settings_txt)
  38 {
  39   QStringList items = settings_txt.split(';');
  40   if (items.size() != 2) {
  41     EG_ERR_RETURN("syntax error for error weighting");
  42   }
  43   m_Err0 = items[0].trimmed().toDouble();
  44   m_ErrS = 0.5*m_Err0;
  45   m_XS   = items[1].trimmed().toDouble();
  46   m_Exp = 1.0;
  47   if (m_Err0 > 1e-3) {
  48     m_Exp = log(m_ErrS/m_Err0)/log(1.0 - m_XS);
  49   }
  50 }
  51
  52 double ErrorFunction::operator()(double x)
  53 {
  54   double e = 0;
  55   if (m_Active) {
  56     double x0 = fabs(1-x);
  57     e = m_Err0*pow(x0, m_Exp);
  58     m_MaxErr = max(m_MaxErr, x0);
  59     m_TotalError += e;
  60     m_AverageError += x0;
  61     ++m_NumCalls;
  62     ++m_NumTotalCalls;
  63   }
  64   return e;
  65 }
  66
  67 double ErrorFunction::averageError()
  68 {
  69   if (m_NumTotalCalls == 0) {
  70     return 0;
  71   }
  72   return m_AverageError/m_NumTotalCalls;
  73 }
  74
  75 double ErrorFunction::totalError()
  76 {
  77   if (m_NumCalls == 0) {
  78     return 0;
  79   }
  80   return m_TotalError/m_NumCalls;
  81 }
  82
  83 void ErrorFunction::reset(bool reset_average)
  84 {
  85   m_MaxErr = 0;
  86   m_TotalError = 0;
  87   m_NumCalls = 0;
  88   if (reset_average) {
  89     m_AverageError = 0;
  90     m_NumTotalCalls = 0;
  91   }
  92 }
  93
  94 Optimisation::Optimisation()
  95 {
  96   setDeltas(1e-6);
  97   F = new double**[3];
  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];
 102     }
 103   }
 104   m_ErrorFunctions.clear();
 105 }
 106
 107 void Optimisation::getErrSet(QString group, QString key, double err0, double xs, ErrorFunction* err_func)
 108 {
 109   QString err0_txt, xs_txt;
 110   err0_txt.setNum(err0);
 111   xs_txt.setNum(xs);
 112   QString value = err0_txt + "; " + xs_txt;
 113   QString variable;
 114   QSettings *qset = GuiMainWindow::settings();
 115   QString typed_key = "string/" + key;
 116   if (group != QObject::tr("General")) {
 117     qset->beginGroup(group);
 118   }
 119   if (!qset->contains(typed_key)) {
 120     qset->setValue(typed_key, value);
 121   }
 122   variable = (qset->value(typed_key,variable)).toString();
 123   if (group != QObject::tr("General")) {
 124     qset->endGroup();
 125   }
 126   err_func->set(variable);
 127   err_func->setName(key.replace(" ", "_"));
 128   m_ErrorFunctions.append(err_func);
 129 }
 130
 131 double Optimisation::angleX(const vec3_t &v1, const vec3_t &v2)
 132 {
 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);
 138 }
 139
 140 void Optimisation::computeDerivatives(vec3_t x)
 141 {
 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);
 169
 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);
 173
 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);
 177
 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);
 181   J[1][0] = J[0][1];
 182   J[2][0] = J[0][2];
 183   J[2][1] = J[1][2];
 184 };
 185
 186 vec3_t Optimisation::optimise(vec3_t x)
 187 {
 188   computeDerivatives(x);
 189   mat3_t M = J.inverse();
 190   return (-1)*(M*grad_f);
 191 };
 192
 193 void Optimisation::resetErrorFunctions(bool reset_average)
 194 {
 195   foreach (ErrorFunction *err_func, m_ErrorFunctions) {
 196     err_func->reset(reset_average);
 197   }
 198 }
 199
 200 double Optimisation::totalError()
 201 {
 202   double e = 0;
 203   foreach (ErrorFunction *err_func, m_ErrorFunctions) {
 204     e += err_func->totalError();
 205   }
 206   return e;
 207 }
 208
 209 void Optimisation::printErrors()
 210 {
 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;
 214     }
 215   }
 216 }