1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
8 //                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
9 //
10 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
11 //////////////////////////////////////////////////////////////////////////////////////
12 
13 
14 #ifndef QMCPLUSPLUS_CUBIC_B_SPLINE_GRID_H
15 #define QMCPLUSPLUS_CUBIC_B_SPLINE_GRID_H
16 #include "Numerics/GridTraits.h"
17 #include "Numerics/BsplineOneDimSolvers.h"
18 #include <limits>
19 
20 /** CubicBsplineGrid
21  *
22  * Empty declaration to be specialized. Three template parameters are
23  * - T data type
24  * - GRIDTYPE enumeration of the grid type
25  * - PBC true for periodic boundary condition
26  */
27 template<class T, unsigned GRIDTYPE, unsigned BC>
28 struct CubicBsplineGrid
29 {};
30 
31 /** specialization for linear grid with PBC */
32 template<class T>
33 struct CubicBsplineGrid<T, LINEAR_1DGRID, PBC_CONSTRAINTS>
34 {
35   typedef typename GridTraits<T>::point_type point_type;
36   typedef typename GridTraits<T>::value_type value_type;
37   typedef std::vector<T> container_type;
38   int i0, i1, i2, i3;
39   point_type GridStart, GridEnd, GridDelta, GridDeltaInv, GridDeltaInv2, L, Linv;
40   point_type curPoint;
41   point_type tp[4];
42 
43   inline CubicBsplineGrid() : curPoint(-10000) {}
44 
45   inline bool getGridPoint(point_type x, int& i)
46   {
47     //point_type delta = x - GridStart;
48     //delta -= std::floor(delta*Linv)*L;
49     //point_type ipart;
50     //point_type t = modf (delta*GridDeltaInv, &ipart);
51     //int i = (int) ipart;
52     point_type delta = x - GridStart;
53     delta -= std::floor(delta * Linv) * L;
54     i            = static_cast<int>(delta * GridDeltaInv);
55     point_type t = delta * GridDeltaInv - i;
56     tp[0]        = t * t * t;
57     tp[1]        = t * t;
58     tp[2]        = t;
59     tp[3]        = 1.0;
60     return true;
61   }
62 
63   void spline(point_type start, point_type end, const container_type& data, container_type& p, bool closed)
64   {
65     GridStart = start;
66     GridEnd   = end;
67     int N     = data.size();
68     if (closed)
69       N--;
70     p.resize(N + 3);
71     L             = end - start;
72     Linv          = 1.0 / L;
73     GridDelta     = L / static_cast<T>(N);
74     GridDeltaInv  = 1.0 / GridDelta;
75     GridDeltaInv2 = 1.0 / GridDelta / GridDelta;
76     SolvePeriodicInterp1D<T>::apply(data, p, N);
77     p[0]     = p[N];
78     p[N + 1] = p[1];
79     p[N + 2] = p[2];
80   }
81 };
82 
83 /** specialization for linear grid with PBC */
84 template<class T>
85 struct CubicBsplineGrid<T, LINEAR_1DGRID, FIRSTDERIV_CONSTRAINTS>
86 {
87   typedef typename GridTraits<T>::point_type point_type;
88   typedef typename GridTraits<T>::value_type value_type;
89   typedef std::vector<T> container_type;
90   typedef std::size_t size_t;
91   ///number of points
92   int Npts;
93   point_type GridStart, GridEnd, GridDelta, GridDeltaInv, GridDeltaInv2, L, Linv;
94   point_type tp[4];
95 
96   inline CubicBsplineGrid() {}
97 
98   inline bool getGridPoint(point_type x, int& i)
99   {
100     if (x < GridStart || x > GridEnd)
101       return false;
102     point_type delta = x - GridStart;
103     delta -= std::floor(delta * Linv) * L;
104     i            = static_cast<int>(delta * GridDeltaInv);
105     point_type t = delta * GridDeltaInv - i;
106     tp[0]        = t * t * t;
107     tp[1]        = t * t;
108     tp[2]        = t;
109     tp[3]        = 1.0;
110     return true;
111   }
112 
113   /** set linear grid
114    * @param start starting grid
115    * @param end ending grid
116    * @param n size of data
117    */
118   void setGrid(point_type start, point_type end, size_t n)
119   {
120     Npts          = n;
121     L             = end - start;
122     Linv          = 1.0 / L;
123     GridDelta     = L / static_cast<point_type>(Npts - 1);
124     GridDeltaInv  = 1.0 / GridDelta;
125     GridDeltaInv2 = 1.0 / GridDelta / GridDelta;
126     GridStart     = start;
127     GridEnd       = end;
128   }
129 
130   void spline(point_type start, point_type end, const container_type& data, container_type& p, bool closed)
131   {
132     setGrid(start, end, data.size());
133     p.resize(Npts + 2);
134     point_type bcLower[] = {-3.0, 0.0, 3.0, 0.0};
135     point_type bcUpper[] = {-3.0, 0.0, 3.0, 0.0};
136     bcLower[3]           = data[1] - data[0];
137     bcUpper[3]           = data[Npts - 1] - data[Npts - 2];
138     SolveFirstDerivInterp1D<point_type>::apply(data, p, Npts, bcLower, bcUpper);
139   }
140 
141   void spline(point_type start,
142               point_type end,
143               value_type startDeriv,
144               value_type endDeriv,
145               const container_type& data,
146               container_type& p)
147   {
148     setGrid(start, end, data.size());
149     p.resize(Npts + 2);
150     point_type bcLower[] = {-3.0, 0.0, 3.0, 0.0};
151     point_type bcUpper[] = {-3.0, 0.0, 3.0, 0.0};
152     bcLower[3]           = startDeriv * GridDelta;
153     bcUpper[3]           = endDeriv * GridDelta;
154     SolveFirstDerivInterp1D<point_type>::apply(data, p, Npts, bcLower, bcUpper);
155   }
156 };
157 #endif
158