Add a SplineFunction class for cubic spline interpolation

Author: robertodr

src/bin/plot_green_spherical.cpp

@@ -14,7 +14,6 @@
 #include "UniformDielectric.hpp"
 
 std::vector<double> generate_grid(double, double, size_t);
-std::string header();
 void case1();
 void case2();
 void case3();
@@ -40,59 +39,48 @@ std::vector<double> generate_grid(double min_val, double max_val, size_t nPoints
     return grid;
 }
 
-std::string header()
-{
-    std::stringstream tmp;
-    tmp << "#" << '\t' << "x_2" << '\t' << "z_2" << '\t';
-    tmp << "gf_value" << '\t' << "image" << '\t' << "Coulomb" << '\t';
-    tmp << "derivativeProbe" << '\t' << "der_image" << '\t' << "der_Coulomb" << std::endl;
-    return tmp.str();
-}
-
 void case1()
 {
-    size_t nPoints = 100;
+    size_t nPoints = 200;
     double epsInside = 80.0;
     double epsOutside = 2.0;
     Eigen::Vector3d sphereCenter = Eigen::Vector3d::Zero();
     double sphereRadius = 100.0;
-    double width = 10.0;
+    double width = 9.0;
 
-    Eigen::Vector3d source;
-    source << 4.0, 0.0, 85.0;
+    Eigen::Vector3d source = Eigen::Vector3d::Zero();
     Eigen::Vector3d probe = Eigen::Vector3d::Zero();
 
-    double x_min = 10.0;
-    double x_max = 300.0;
-    std::vector<double> x = generate_grid(x_min, x_max, nPoints);
-    double z_min = 10.0;
-    double z_max = 300.0;
+    double z_min = -20.0;
+    double z_max = 20.0;
     std::vector<double> z = generate_grid(z_min, z_max, nPoints);
+    double delta = 0.1;
 
     std::ofstream out;
 
     SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter);
     LOG(gf);
 
     out.open("gf_spherical_CASE1.dat");
-    out << header();
     out.precision(16);
     for (size_t i = 0; i < nPoints; ++i) {
+	source << 0.0, 0.0, z[i];
         for (size_t j = 0; j < nPoints; ++j) {
-            probe << x[i], 0.0, z[j];
-            out << '\t' << x[i]
+            probe << 0.0, 0.0, (z[j] + delta);
+            out << '\t' << i
+		<< '\t' << j
+		<< '\t' << z[i]
                 << '\t' << z[j]
-                << '\t' << gf.kernelS(source, probe)
-                << '\t' << gf.imagePotential(source, probe)
                 << '\t' << gf.Coulomb(source, probe)
-                << '\t' << gf.derivativeProbe(Eigen::Vector3d::UnitZ(), source, probe)
-                << '\t' << gf.CoulombDerivative(Eigen::Vector3d::UnitZ(), source, probe)
-                << '\t' << gf.imagePotentialDerivative(Eigen::Vector3d::UnitZ(), source, probe) << std::endl;
+                << '\t' << gf.imagePotential(source, probe)
+                << '\t' << gf.kernelS(source, probe)
+		<< '\t' << (1.0 / gf.coefficientCoulomb(source, probe)) << std::endl;
         }
     }
     out.close();
     LOG_TIME;
 
+    /*
     UniformDielectric<AD_directional, CollocationIntegrator> gf_inside(epsInside);
     out.open("gf_uniform_inside_CASE1.dat");
     out << "#" << '\t' << "x_2" << '\t' << "z_2" << '\t';
@@ -124,4 +112,5 @@ void case1()
         }
     }
     out.close();
+    */
 }

src/utils/SplineFunction.hpp

@@ -0,0 +1,78 @@
+/* pcmsolver_copyright_start */
+/*
+ *     PCMSolver, an API for the Polarizable Continuum Model
+ *     Copyright (C) 2013 Roberto Di Remigio, Luca Frediani and contributors
+ *
+ *     This file is part of PCMSolver.
+ *
+ *     PCMSolver is free software: you can redistribute it and/or modify
+ *     it under the terms of the GNU Lesser General Public License as published by
+ *     the Free Software Foundation, either version 3 of the License, or
+ *     (at your option) any later version.
+ *
+ *     PCMSolver is distributed in the hope that it will be useful,
+ *     but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *     GNU Lesser General Public License for more details.
+ *
+ *     You should have received a copy of the GNU Lesser General Public License
+ *     along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *     For information on the complete list of contributors to the
+ *     PCMSolver API, see: <http://pcmsolver.github.io/pcmsolver-doc>
+ */
+/* pcmsolver_copyright_end */
+
+#ifndef SPLINEFUNCTION_HPP
+#define SPLINEFUNCTION_HPP
+
+#include "Config.hpp"
+
+#include <Eigen/Core>
+#include <unsupported/Eigen/Splines>
+
+/*! \file SplineFunction.hpp
+ *  \class SplineFunction
+ *  \brief Spline interpolation of a function
+ *  \author Roberto Di Remigio
+ *  \date 2015
+ *
+ *  Taken from StackOverflow http://stackoverflow.com/a/29825204/2528668
+ */
+   
+/*! Scale abscissa value to [0, 1] */
+inline double scale(double x) {
+	return (x - xMin_) / (xMax_ - xMin_);
+}
+
+/*! Scale abscissa values to [0, 1] */
+inline Eigen::RowVectorXd scale(const Eigen::VectorXd & x_vec) {
+    return x_vec.unaryExpr([](double x) -> double { return scale(x); }).transpose();
+}
+
+class SplineFunction
+{
+private:
+    typedef Eigen::Spline<double, 1, 3> SplineFit;
+public:
+    /*! \param[in] x vector with abscissa values 
+     *  \param[in] y vector with function values
+     */
+    SplineFunction(const Eigen::VectorXd & x, const Eigen::VectorXd & y) 
+	    : xMin_(x.minCoeff()), xMax_(x.maxCoeff()),
+	    spline_(Eigen::SplineFitting<SplineFit>::Interpolate(y.transpose(), 
+		      std::min<int>(x.rows() - 1, 3), scale(x))
+	           )
+	{}
+    /*! \brief Evaluate spline at given point
+     *  \param[in] x evaluation point
+     */
+    double operator()(double x) const {
+	return spline_(scale(x))(0);
+    }
+private:
+    double xMin_;
+    double xMax_;
+    SplineFit spline_;
+};
+#endif // SPLINEFUNCTION_HPP