Fixed scaling problem by preconditioning with scales based on parametrization gradients

contrib/HighOrderMeshOptimizer/OptHOM.cpp

@@ -149,45 +149,82 @@ void printProgressFunc(const alglib::real_1d_array &x, double Obj, void *HOInst)
 
 
 
+void OptHOM::calcScale(alglib::real_1d_array &scale)
+{
+
+  scale.setlength(mesh.nPC());
+
+  // Calculate scale
+  for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
+    std::vector<double> scaleFV(mesh.nPCFV(iFV),1.);
+    mesh.pcScale(iFV,scaleFV);
+    for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++) scale[mesh.indPCFV(iFV,iPC)] = scaleFV[iPC];
+  }
+
+  // Normalize scale vector (otherwise ALGLIB routines may fail)
+  double scaleNormSq = 0.;
+  for (int i = 0; i < mesh.nPC(); i++) scaleNormSq += scale[i]*scale[i];
+  const double scaleNorm = sqrt(scaleNormSq);
+  for (int i = 0; i < mesh.nPC(); i++) scale[i] /= scaleNorm;
+
+}
+
+
+
+
 void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &initGradObj, int itMax)
 {
 
   static const double EPSG = 0.;
   static const double EPSF = 0.;
   static const double EPSX = 0.;
+  static int OPTMETHOD = 1;
 
   std::cout << "--- Optimization pass with jacBar = " << jacBar << ", lambda = " << lambda << ", lambda2 = " << lambda2 << std::endl;
 
-//  alglib::minlbfgsstate state;
-//  alglib::minlbfgsreport rep;
-  alglib::mincgstate state;
-  alglib::mincgreport rep;
-
-//  minlbfgscreate(3, x, state);
-//  minlbfgssetcond(state, EPSG, EPSF, EPSX, itMax);
-//  minlbfgssetxrep(state, true);
-  mincgcreate(x, state);
-  mincgsetcond(state, EPSG, EPSF, EPSX, itMax);
-  mincgsetxrep(state, true);
-
   iter = 0;
 
-//  alglib::minlbfgsoptimize(state, evalObjGradFunc, printProgressFunc, this);
-  alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this);
-
-//  minlbfgsresults(state, x, rep);
-  mincgresults(state, x, rep);
+  int iterationscount = 0, nfev = 0, terminationtype = -1;
+  if (OPTMETHOD == 1) {
+    alglib::mincgstate state;
+    alglib::mincgreport rep;
+    mincgcreate(x, state);
+    alglib::real_1d_array scale;
+    calcScale(scale);
+    mincgsetscale(state,scale);
+    mincgsetprecscale(state);
+    mincgsetcond(state, EPSG, EPSF, EPSX, itMax);
+    mincgsetxrep(state, true);
+    alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this);
+    mincgresults(state, x, rep);
+    iterationscount = rep.iterationscount;
+    nfev = rep.nfev;
+    terminationtype = rep.terminationtype;
+  }
+  else {
+    alglib::minlbfgsstate state;
+    alglib::minlbfgsreport rep;
+    minlbfgscreate(3, x, state);
+    alglib::real_1d_array scale;
+    calcScale(scale);
+    minlbfgssetscale(state,scale);
+    minlbfgssetprecscale(state);
+    minlbfgssetcond(state, EPSG, EPSF, EPSX, itMax);
+    minlbfgssetxrep(state, true);
+    alglib::minlbfgsoptimize(state, evalObjGradFunc, printProgressFunc, this);
+    minlbfgsresults(state, x, rep);
+    iterationscount = rep.iterationscount;
+    nfev = rep.nfev;
+    terminationtype = rep.terminationtype;
+  }
 
-  std::cout << "Optimization finalized after " << rep.iterationscount << " iterations (" << rep.nfev << " functions evaluations)";
-  switch(int(rep.terminationtype)) {
-//  case -2: std::cout << ", because rounding errors prevented further improvement"; break;
-//  case -1: std::cout << ", because incorrect parameters were specified"; break;
+  std::cout << "Optimization finalized after " << iterationscount << " iterations (" << nfev << " functions evaluations)";
+  switch(int(terminationtype)) {
   case 1: std::cout << ", because relative function improvement is no more than EpsF"; break;
   case 2: std::cout << ", because relative step is no more than EpsX"; break;
   case 4: std::cout << ", because gradient norm is no more than EpsG"; break;
   case 5: std::cout << ", because the maximum number of steps was taken"; break;
-//  case 7: std::cout << ", because stopping conditions are too stringent, further improvement is impossible"; break;
-  default: std::cout << " with code " << int(rep.terminationtype); break;
+  default: std::cout << " with code " << int(terminationtype); break;
   }
   std::cout << "." << std::endl;
 

contrib/HighOrderMeshOptimizer/OptHOM.h

@@ -39,11 +39,12 @@ private:
   double initObj, initMaxDist, initAvgDist;  // Values for reporting
   double minJac, maxJac, maxDist, avgDist;  // Values for reporting
 
-  inline double setBarrierTerm(double jacBarrier) { bTerm = jacBarrier/(1.-jacBarrier); }
+  inline void setBarrierTerm(double jacBarrier) { bTerm = jacBarrier/(1.-jacBarrier); }
   inline double compute_f(double v);
   inline double compute_f1(double v);
   bool addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj);
   bool addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real_1d_array &gradObj);
+  void calcScale(alglib::real_1d_array &scale);
   void OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &initGradObj, int itMax);
 
 };
@@ -83,7 +84,6 @@ inline double OptHOM::compute_f1(double v)
 void OptHOM::getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD)
 {
   minJ = minJac; maxJ = maxJac; maxD = maxDist; avgD = avgDist;
-  std::cout << "getJacDist: minJac = " << minJac << ", maxJac = " << maxJac << ", maxDist = " << maxDist << ", avgDist = " << avgDist << std::endl;
 }
 
 

contrib/HighOrderMeshOptimizer/OptHomMesh.cpp

@@ -440,6 +440,25 @@ void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ)
 
 
 
+void Mesh::pcScale(int iFV, std::vector<double> &scale)
+{
+
+  // Calc. derivative of x, y & z w.r.t. parametric coordinates
+  const SPoint3 dX(1.,0.,0.), dY(0.,1.,0.), dZ(0.,0.,1.);
+  SPoint3 gX, gY, gZ;
+  _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dX,gX);
+  _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dY,gY);
+  _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dZ,gZ);
+
+  // Scale = inverse norm. of vector (dx/du, dy/du, dz/du)
+  scale[0] = 1./sqrt(gX[0]*gX[0]+gY[0]*gY[0]+gZ[0]*gZ[0]);
+  if (_nPCFV[iFV] >= 2) scale[1] = 1./sqrt(gX[1]*gX[1]+gY[1]*gY[1]+gZ[1]*gZ[1]);
+  if (_nPCFV[iFV] == 3) scale[2] = 1./sqrt(gX[2]*gX[2]+gY[2]*gY[2]+gZ[2]*gZ[2]);
+
+}
+
+
+
 void Mesh::writeMSH(const char *filename)
 {
   FILE *f = fopen(filename, "w");

contrib/HighOrderMeshOptimizer/OptHomMesh.h

@@ -37,6 +37,8 @@ public:
   inline double distSq(int iFV);
   inline void gradDistSq(int iV, std::vector<double> &gDSq);
 
+  void pcScale(int iFV, std::vector<double> &scale);          // Calc. scale of parametric coordinates for preconditioning
+
   void getUvw(double *it);
   void updateMesh(const double *it);
   void distSqToStraight(std::vector<double> &dSq);