1 /*
2  * Copyright 2010-2020 The VOTCA Development Team (http://www.votca.org)
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  *     http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  *
16  */
17 
18 #ifndef VOTCA_TOOLS_SPLINE_H
19 #define VOTCA_TOOLS_SPLINE_H
20 
21 // Local VOTCA includes
22 #include "eigen.h"
23 #include "types.h"
24 
25 namespace votca {
26 namespace tools {
27 
28 /**
29  * \brief Spline Class
30  *
31  *  class supports spline interpolation and fit of data
32  *  the cubic spline class, akima spline class and linear spline class are
33  * inherited from this one
34  */
35 
36 class Spline {
37  public:
38   Spline() = default;
39 
40   virtual ~Spline() = default;
41 
42   /**
43    * \brief Calculate interpolating spline for given (x,y) values. Points on
44    * resulting spline can be obtained via Calculate(). \param x values of data
45    * to be interpolated \param y values of data to be interpolated both vectors
46    * must be of same size
47    */
48   virtual void Interpolate(const Eigen::VectorXd &x,
49                            const Eigen::VectorXd &y) = 0;
50 
51   /**
52    * \brief Fit spline through noisy (x,y) values. Points on resulting fitted
53    * spline can be obtained via Calculate(). \param x values of data to be
54    * fitted \param y values of data to be fitted both vectors must be of same
55    * size
56    */
57   virtual void Fit(const Eigen::VectorXd &x, const Eigen::VectorXd &y) = 0;
58 
59   /**
60    * \brief Calculate spline function value for a given x value on the spline
61    * created by Interpolate() or Fit() \param x data value \return y value
62    */
63   virtual double Calculate(double x) = 0;
64 
65   /**
66    * \brief Calculate y value for a given x value on the derivative of the
67    * spline created by function Interpolate or Fit \param x data value \return y
68    * value of derivative
69    */
70   virtual double CalculateDerivative(double x) = 0;
71 
72   /// enum for type of boundary condition
73   enum eBoundary {
74     splineNormal = 0,         ///< normal boundary conditions: \f$f_0=f_N=0\f$
75     splinePeriodic = 1,       ///< periodic boundary conditions: \f$f_0=f_N\f$
76     splineDerivativeZero = 2  ///< derivatives and end-points are zero.
77   };
78 
79   /**
80    * \brief Set the boundary type of the spline
81    * \param boundary of type eBoundary
82    */
83   void setBC(eBoundary bc) { boundaries_ = bc; }
84 
85   /**
86    * \brief Set the boundary type of the spline
87    * \param boundary of type int
88    */
89   void setBCInt(Index bc) {
90     switch (bc) {
91       case 0:
92         boundaries_ = splineNormal;
93         break;
94       case 1:
95         boundaries_ = splinePeriodic;
96         break;
97       case 2:
98         boundaries_ = splineDerivativeZero;
99         break;
100     }
101   }
102 
103   /**
104    * \brief Get the grid point of certain index
105    * \param index of grid point
106    * \return grid value
107    */
108   double getGridPoint(int i);
109 
110   /**
111    * \brief Calculate spline function values for given x values on the spline
112    * created by Interpolate() or Fit() \param vector of x data values \return
113    * vector of y value
114    */
115   Eigen::VectorXd Calculate(const Eigen::VectorXd &x);
116 
117   /**
118    * \brief Calculate y values for given x values on the derivative of the
119    * spline created by function Interpolate or Fit \param vector of x data
120    * values \return vector of y value
121    */
122   Eigen::VectorXd CalculateDerivative(const Eigen::VectorXd &x);
123 
124   /**
125    * \brief Print spline values (using Calculate()) on output "out" on the
126    * entire grid in steps of "interval" \param reference "out" to output \param
127    * steps of size "interval"
128    */
129   void Print(std::ostream &out, double interval);
130 
131   /**
132    * \brief Determine the index of the interval containing value r
133    * \param value r
134    * \return interval index
135    */
136   Index getInterval(double r);
137 
138   /**
139    * \brief Generate the grid for fitting from "min" to "max" in steps of "h"
140    * \param left interval border "min"
141    * \param right interval border "max"
142    * \param step "h"
143    * \return number of grid values in the interval
144    */
145   Index GenerateGrid(double min, double max, double h);
146 
147   /**
148    * \brief Get the grid array x
149    * \return pointer to the corresponding array
150    */
151   Eigen::VectorXd &getX() { return r_; }
152   const Eigen::VectorXd &getX() const { return r_; }
153   /**
154    * \brief Get the spline data  f_
155    * \return reference to the corresponding array
156    */
157   // Eigen::VectorXd &getSplineF() { return  f_; }
158   // const Eigen::VectorXd &getSplineF() const { return  f_; }
159 
160   /**
161    * \brief Get second derivatives (cubic splines)
162    * \return reference to the corresponding array
163    */
164   // Eigen::VectorXd &getSplineF2() { return  f2_; }
165   // const Eigen::VectorXd &getSplineF2() const { return  f_; }
166 
167  protected:
168   eBoundary boundaries_ = eBoundary::splineNormal;
169   // the grid points
170   Eigen::VectorXd r_;
171 };
172 
173 }  // namespace tools
174 }  // namespace votca
175 
176 #endif  // VOTCA_TOOLS_SPLINE_H
177