some simple global optimization routines

Former-commit-id: f12871ad8ef24396a8d8431e1ff0b8c8c34f6653

include/igl/grid_search.cpp

@@ -0,0 +1,64 @@
+#include "grid_search.h"
+#include <iostream>
+#include <cassert>
+
+template <
+  typename Scalar, 
+  typename DerivedX, 
+  typename DerivedLB, 
+  typename DerivedUB, 
+  typename DerivedI>
+IGL_INLINE Scalar igl::grid_search(
+  const std::function< Scalar (DerivedX &) > f,
+  const Eigen::MatrixBase<DerivedLB> & LB,
+  const Eigen::MatrixBase<DerivedUB> & UB,
+  const Eigen::MatrixBase<DerivedI> & I,
+  DerivedX & X)
+{
+  Scalar fval = std::numeric_limits<Scalar>::max();
+  const int dim = LB.size();
+  assert(UB.size() == dim && "UB should match LB size");
+  assert(I.size() == dim && "I should match LB size");
+  X.resize(dim);
+
+  // Working X value
+  DerivedX Xrun(dim);
+  std::function<void(const int, DerivedX &)> looper;
+  int calls = 0;
+  looper = [&](
+    const int d,
+    DerivedX & Xrun)
+  {
+    assert(d < dim);
+    Eigen::Matrix<Scalar,Eigen::Dynamic,1> vals = 
+      Eigen::Matrix<Scalar,Eigen::Dynamic,1>::LinSpaced(I(d),LB(d),UB(d));
+    for(int c = 0;c<I(d);c++)
+    {
+      Xrun(d) = vals(c);
+      if(d+1 < dim)
+      {
+        looper(d+1,Xrun);
+      }else
+      {
+        //std::cout<<"call: "<<calls<<std::endl;
+        // Base case
+        const Scalar val = f(Xrun);
+        calls++;
+        if(val < fval)
+        {
+          fval = val;
+          X = Xrun;
+          std::cout<<calls<<": "<<fval<<" | "<<X<<std::endl;
+        }
+      }
+    }
+  };
+  looper(0,Xrun);
+
+  return fval;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template double igl::grid_search<double, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<int, 1, 3, 1, 1, 3> >(std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3>&)>, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, 1, 3, 1, 1, 3> > const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>&);
+template float igl::grid_search<float, Eigen::Matrix<float, 1, -1, 1, 1, -1>, Eigen::Matrix<float, 1, -1, 1, 1, -1>, Eigen::Matrix<float, 1, -1, 1, 1, -1>, Eigen::Matrix<int, 1, -1, 1, 1, -1> >(std::function<float (Eigen::Matrix<float, 1, -1, 1, 1, -1>&)>, Eigen::MatrixBase<Eigen::Matrix<float, 1, -1, 1, 1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, -1, 1, 1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, 1, -1, 1, 1, -1> > const&, Eigen::Matrix<float, 1, -1, 1, 1, -1>&);
+#endif

include/igl/grid_search.h

@@ -0,0 +1,42 @@
+#ifndef IGL_GRID_SEARCH_H
+#define IGL_GRID_SEARCH_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <functional>
+namespace igl
+{
+  // Solve the problem:
+  //
+  //   minimize f(x)
+  //   subject to lb ≤ x ≤ ub 
+  // 
+  // by exhaustive grid search.
+  //
+  // Inputs:
+  //   f  function to minimize
+  //   LB  #X vector of finite lower bounds
+  //   UB  #X vector of finite upper bounds
+  //   I  #X vector of number of steps for each variable
+  // Outputs:
+  //   X  #X optimal parameter vector
+  // Returns f(X)
+  //
+  template <
+    typename Scalar, 
+    typename DerivedX, 
+    typename DerivedLB, 
+    typename DerivedUB, 
+    typename DerivedI>
+  IGL_INLINE Scalar grid_search(
+    const std::function< Scalar (DerivedX &) > f,
+    const Eigen::MatrixBase<DerivedLB> & LB,
+    const Eigen::MatrixBase<DerivedUB> & UB,
+    const Eigen::MatrixBase<DerivedI> & I,
+    DerivedX & X);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "grid_search.cpp"
+#endif
+
+#endif

include/igl/pso.cpp

@@ -0,0 +1,174 @@
+#include "pso.h"
+#include <cassert>
+#include <Eigen/StdVector>
+#include <vector>
+#include <iostream>
+
+template <
+  typename Scalar, 
+  typename DerivedX,
+  typename DerivedLB, 
+  typename DerivedUB>
+IGL_INLINE Scalar igl::pso(
+  const std::function< Scalar (DerivedX &) > f,
+  const Eigen::MatrixBase<DerivedLB> & LB,
+  const Eigen::MatrixBase<DerivedUB> & UB,
+  const int max_iters,
+  const int population,
+  DerivedX & X)
+{
+  const Eigen::Array<bool,Eigen::Dynamic,1> P =
+    Eigen::Array<bool,Eigen::Dynamic,1>::Zero(LB.size(),1);
+  return igl::pso(f,LB,UB,P,max_iters,population,X);
+}
+
+template <
+  typename Scalar, 
+  typename DerivedX,
+  typename DerivedLB, 
+  typename DerivedUB,
+  typename DerivedP>
+IGL_INLINE Scalar igl::pso(
+  const std::function< Scalar (DerivedX &) > f,
+  const Eigen::MatrixBase<DerivedLB> & LB,
+  const Eigen::MatrixBase<DerivedUB> & UB,
+  const Eigen::DenseBase<DerivedP> & P,
+  const int max_iters,
+  const int population,
+  DerivedX & X)
+{
+  const int dim = LB.size();
+  assert(UB.size() == dim && "UB should match LB size");
+  assert(P.size() == dim && "P should match LB size");
+  typedef std::vector<DerivedX,Eigen::aligned_allocator<DerivedX> > VectorList;
+  VectorList position(population);
+  VectorList best_position(population);
+  VectorList velocity(population);
+  Eigen::Matrix<Scalar,Eigen::Dynamic,1> best_f(population);
+  // https://en.wikipedia.org/wiki/Particle_swarm_optimization#Algorithm
+  //
+  // g → X
+  // p_i → best[i]
+  // v_i → velocity[i]
+  // x_i → position[i]
+  Scalar min_f = std::numeric_limits<Scalar>::max();
+  for(int p=0;p<population;p++)
+  {
+    {
+      const DerivedX R = DerivedX::Random(dim).array()*0.5+0.5;
+      position[p] = LB.array() + R.array()*(UB-LB).array();
+    }
+    best_f[p] = f(position[p]);
+    best_position[p] = position[p];
+    if(best_f[p] < min_f)
+    {
+      min_f = best_f[p];
+      X = best_position[p];
+    }
+    {
+      const DerivedX R = DerivedX::Random(dim);
+      velocity[p] = (UB-LB).array() * R.array();
+    }
+  }
+
+  int iter = 0;
+  Scalar omega = 0.98;
+  Scalar phi_p = 0.01;
+  Scalar phi_g = 0.01;
+  while(true)
+  {
+    //if(iter % 10 == 0)
+    //{
+    //  std::cout<<iter<<":"<<std::endl;
+    //  for(int p=0;p<population;p++)
+    //  {
+    //    std::cout<<"  "<<best_f[p]<<", "<<best_position[p]<<std::endl;
+    //  }
+    //  std::cout<<std::endl;
+    //}
+
+    for(int p=0;p<population;p++)
+    {
+      const DerivedX R_p = DerivedX::Random(dim).array()*0.5+0.5;
+      const DerivedX R_g = DerivedX::Random(dim).array()*0.5+0.5;
+      velocity[p] = 
+        omega * velocity[p].array() +
+        phi_p * R_p.array() *(best_position[p] - position[p]).array() + 
+        phi_g * R_g.array() *(               X - position[p]).array();
+      position[p] += velocity[p];
+      // Clamp to bounds
+      for(int d = 0;d<dim;d++)
+      {
+//#define IGL_PSO_REFLECTION
+#ifdef IGL_PSO_REFLECTION
+        assert(!P(d));
+        // Reflect velocities if exceeding bounds
+        if(position[p](d) < LB(d))
+        {
+          position[p](d) = LB(d);
+          if(velocity[p](d) < 0.0) velocity[p](d) *= -1.0;
+        }
+        if(position[p](d) > UB(d))
+        {
+          position[p](d) = UB(d);
+          if(velocity[p](d) > 0.0) velocity[p](d) *= -1.0;
+        }
+#else
+//#warning "trying no bounds on periodic"
+//        // TODO: I'm not sure this is the right thing to do/enough. The
+//        // velocities could be weird. Suppose the current "best" value is ε and
+//        // the value is -ε and the "periodic bounds" [0,2π]. Moding will send
+//        // the value to 2π-ε but the "velocity" term will now be huge pointing
+//        // all the way from 2π-ε to ε.
+//        //
+//        // Q: Would it be enough to try (all combinations) of ±(UB-LB) before
+//        // computing velocities to "best"s? In the example above, instead of
+//        //
+//        //     v += best - p = ε - (2π-ε) = -2π+2ε
+//        //
+//        // you'd use
+//        //
+//        //     v +=  / argmin  |b - p|            \  - p = (ε+2π)-(2π-ε) = 2ε
+//        //          |                              |
+//        //           \ b∈{best, best+2π, best-2π} /
+//        //
+//        // Though, for multivariate b,p,v this would seem to explode
+//        // combinatorially.
+//        //
+//        // Maybe periodic things just shouldn't be bounded and we hope that the
+//        // forces toward the current minima "regularize" them away from insane
+//        // values.
+//        if(P(d))
+//        {
+//          position[p](d) = std::fmod(position[p](d)-LB(d),UB(d)-LB(d))+LB(d);
+//        }else
+//        {
+//          position[p](d) = std::max(LB(d),std::min(UB(d),position[p](d)));
+//        }
+        position[p](d) = std::max(LB(d),std::min(UB(d),position[p](d)));
+#endif
+      }
+      const Scalar fp = f(position[p]);
+      if(fp<best_f[p])
+      {
+        best_f[p] = fp;
+        best_position[p] = position[p];
+        if(best_f[p] < min_f)
+        {
+          min_f = best_f[p];
+          X = best_position[p];
+        }
+      }
+    }
+    iter++;
+    if(iter>=max_iters)
+    {
+      break;
+    }
+  }
+  return min_f;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template float igl::pso<float, Eigen::Matrix<float, 1, -1, 1, 1, -1>, Eigen::Matrix<float, 1, -1, 1, 1, -1>, Eigen::Matrix<float, 1, -1, 1, 1, -1> >(std::function<float (Eigen::Matrix<float, 1, -1, 1, 1, -1>&)>, Eigen::MatrixBase<Eigen::Matrix<float, 1, -1, 1, 1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, -1, 1, 1, -1> > const&, int, int, Eigen::Matrix<float, 1, -1, 1, 1, -1>&);
+#endif

include/igl/pso.h

@@ -0,0 +1,59 @@
+#ifndef IGL_PSO_H
+#define IGL_PSO_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include <functional>
+
+namespace igl
+{
+  // Solve the problem:
+  //
+  //   minimize f(x)
+  //   subject to lb ≤ x ≤ ub 
+  // 
+  // by particle swarm optimization (PSO).
+  //
+  // Inputs:
+  //   f  function that evaluates the objective for a given "particle" location
+  //   LB  #X vector of lower bounds 
+  //   UB  #X vector of upper bounds 
+  //   max_iters  maximum number of iterations
+  //   population  number of particles in swarm
+  // Outputs:
+  //   X  best particle seen so far
+  // Returns objective corresponding to best particle seen so far
+  template <
+    typename Scalar, 
+    typename DerivedX,
+    typename DerivedLB, 
+    typename DerivedUB>
+  IGL_INLINE Scalar pso(
+    const std::function< Scalar (DerivedX &) > f,
+    const Eigen::MatrixBase<DerivedLB> & LB,
+    const Eigen::MatrixBase<DerivedUB> & UB,
+    const int max_iters,
+    const int population,
+    DerivedX & X);
+  // Inputs:
+  //   P  whether each DOF is periodic
+  template <
+    typename Scalar, 
+    typename DerivedX,
+    typename DerivedLB, 
+    typename DerivedUB,
+    typename DerivedP>
+  IGL_INLINE Scalar pso(
+    const std::function< Scalar (DerivedX &) > f,
+    const Eigen::MatrixBase<DerivedLB> & LB,
+    const Eigen::MatrixBase<DerivedUB> & UB,
+    const Eigen::DenseBase<DerivedP> & P,
+    const int max_iters,
+    const int population,
+    DerivedX & X);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "pso.cpp"
+#endif
+
+#endif

include/igl/random_search.cpp

@@ -0,0 +1,36 @@
+#include "random_search.h"
+#include <iostream>
+#include <cassert>
+
+template <
+  typename Scalar, 
+  typename DerivedX, 
+  typename DerivedLB, 
+  typename DerivedUB>
+IGL_INLINE Scalar igl::random_search(
+  const std::function< Scalar (DerivedX &) > f,
+  const Eigen::MatrixBase<DerivedLB> & LB,
+  const Eigen::MatrixBase<DerivedUB> & UB,
+  const int iters,
+  DerivedX & X)
+{
+  Scalar min_f = std::numeric_limits<Scalar>::max();
+  const int dim = LB.size();
+  assert(UB.size() == dim && "UB should match LB size");
+  for(int iter = 0;iter<iters;iter++)
+  {
+    const DerivedX R = DerivedX::Random(dim).array()*0.5+0.5;
+    DerivedX Xr = LB.array() + R.array()*(UB-LB).array();
+    const Scalar fr = f(Xr);
+    if(fr<min_f)
+    {
+      min_f = fr;
+      X = Xr;
+    }
+  }
+  return min_f;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template float igl::random_search<float, Eigen::Matrix<float, 1, -1, 1, 1, -1>, Eigen::Matrix<float, 1, -1, 1, 1, -1>, Eigen::Matrix<float, 1, -1, 1, 1, -1> >(std::function<float (Eigen::Matrix<float, 1, -1, 1, 1, -1>&)>, Eigen::MatrixBase<Eigen::Matrix<float, 1, -1, 1, 1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, -1, 1, 1, -1> > const&, int, Eigen::Matrix<float, 1, -1, 1, 1, -1>&);
+#endif

include/igl/random_search.h

@@ -0,0 +1,42 @@
+#ifndef IGL_RANDOM_SEARCH_H
+#define IGL_RANDOM_SEARCH_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <functional>
+namespace igl
+{
+  // Solve the problem:
+  //
+  //   minimize f(x)
+  //   subject to lb ≤ x ≤ ub 
+  // 
+  // by uniform random search.
+  //
+  // Inputs:
+  //   f  function to minimize
+  //   LB  #X vector of finite lower bounds
+  //   UB  #X vector of finite upper bounds
+  //   iters  number of iterations
+  // Outputs:
+  //   X  #X optimal parameter vector
+  // Returns f(X)
+  //
+  template <
+    typename Scalar, 
+    typename DerivedX, 
+    typename DerivedLB, 
+    typename DerivedUB>
+  IGL_INLINE Scalar random_search(
+    const std::function< Scalar (DerivedX &) > f,
+    const Eigen::MatrixBase<DerivedLB> & LB,
+    const Eigen::MatrixBase<DerivedUB> & UB,
+    const int iters,
+    DerivedX & X);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "random_search.cpp"
+#endif
+
+#endif
+