Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2021 The VES code team
       3             :    (see the PEOPLE-VES file at the root of this folder for a list of names)
       4             : 
       5             :    See http://www.ves-code.org for more information.
       6             : 
       7             :    This file is part of VES code module.
       8             : 
       9             :    The VES code module is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    The VES code module is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with the VES code module.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : 
      23             : #include "GridLinearInterpolation.h"
      24             : 
      25             : #include "tools/Grid.h"
      26             : #include "tools/Exception.h"
      27             : #include "tools/Tools.h"
      28             : 
      29             : 
      30             : namespace PLMD {
      31             : namespace ves {
      32             : 
      33             : 
      34        1869 : double GridLinearInterpolation::getGridValueWithLinearInterpolation_1D(GridBase* grid_pntr, const std::vector<double>& arg) {
      35             : 
      36        1869 :   plumed_massert(grid_pntr->getDimension()==1,"The grid is of the wrong dimension, should be one-dimensional");
      37        1869 :   plumed_massert(arg.size()==1,"input value is of the wrong size");
      38             : 
      39        1869 :   double x = arg[0];
      40        1869 :   double grid_dx = grid_pntr->getDx()[0];
      41        1869 :   double grid_min; Tools::convert( grid_pntr->getMin()[0], grid_min);
      42        1869 :   std::vector<unsigned int> i0(1); i0[0] = unsigned( std::floor( (x-grid_min)/grid_dx ) );
      43        1869 :   std::vector<unsigned int> i1(1); i1[0] = unsigned( std::ceil(  (x-grid_min)/grid_dx ) );
      44             :   //
      45        1869 :   double x0 = grid_pntr->getPoint(i0)[0];
      46        1869 :   double x1 = grid_pntr->getPoint(i1)[0];
      47        1869 :   double y0 = grid_pntr->getValue(i0);
      48        1869 :   double y1 = grid_pntr->getValue(i1);
      49             :   //
      50        1869 :   return linearInterpolation(x,x0,x1,y0,y1);
      51             : }
      52             : 
      53             : 
      54       47164 : double GridLinearInterpolation::getGridValueWithLinearInterpolation_2D(GridBase* grid_pntr, const std::vector<double>& arg) {
      55       47164 :   plumed_massert(grid_pntr->getDimension()==2,"The grid is of the wrong dimension, should be two-dimensional");
      56       47164 :   plumed_massert(arg.size()==2,"input value is of the wrong size");
      57             : 
      58       47164 :   std::vector<double> grid_delta = grid_pntr->getDx();
      59       47164 :   std::vector<double> grid_min(2);
      60       47164 :   Tools::convert( grid_pntr->getMin()[0], grid_min[0]);
      61       47164 :   Tools::convert( grid_pntr->getMin()[1], grid_min[1]);
      62             : 
      63       47164 :   std::vector<unsigned int> i00(2);
      64       47164 :   std::vector<unsigned int> i01(2);
      65       47164 :   std::vector<unsigned int> i10(2);
      66       47164 :   std::vector<unsigned int> i11(2);
      67             : 
      68       47164 :   i00[0] = i01[0] = unsigned( std::floor( (arg[0]-grid_min[0])/grid_delta[0] ) );
      69       47164 :   i10[0] = i11[0] = unsigned( std::ceil(  (arg[0]-grid_min[0])/grid_delta[0] ) );
      70             : 
      71       47164 :   i00[1] = i10[1] = unsigned( std::floor( (arg[1]-grid_min[1])/grid_delta[1] ) );
      72       47164 :   i01[1] = i11[1] = unsigned( std::ceil(  (arg[1]-grid_min[1])/grid_delta[1] ) );
      73             : 
      74             :   // https://en.wikipedia.org/wiki/Bilinear_interpolation
      75       47164 :   double x = arg[0];
      76       47164 :   double y = arg[1];
      77             : 
      78       47164 :   double x0 = grid_pntr->getPoint(i00)[0];
      79       47164 :   double x1 = grid_pntr->getPoint(i10)[0];
      80       47164 :   double y0 = grid_pntr->getPoint(i00)[1];
      81       47164 :   double y1 = grid_pntr->getPoint(i11)[1];
      82             : 
      83       47164 :   double f00 = grid_pntr->getValue(i00);
      84       47164 :   double f01 = grid_pntr->getValue(i01);
      85       47164 :   double f10 = grid_pntr->getValue(i10);
      86       47164 :   double f11 = grid_pntr->getValue(i11);
      87             : 
      88             :   // linear interpolation in x-direction
      89             :   double fx0 = linearInterpolation(x,x0,x1,f00,f10);
      90             :   double fx1 = linearInterpolation(x,x0,x1,f01,f11);
      91             :   // linear interpolation in y-direction
      92             :   double fxy = linearInterpolation(y,y0,y1,fx0,fx1);
      93             :   //
      94       47164 :   return fxy;
      95             : }
      96             : 
      97             : 
      98           0 : double GridLinearInterpolation::getGridValueWithLinearInterpolation_ND(GridBase* grid_pntr, const std::vector<double>& arg) {
      99           0 :   unsigned int dimension = grid_pntr->getDimension();
     100           0 :   plumed_massert(dimension==arg.size(),"The grid dimensions do not match the the given arguments.");
     101             :   // get points first
     102           0 :   std::vector<std::vector<unsigned>> point_indices = GridLinearInterpolation::getAdjacentPoints(grid_pntr, arg);
     103             : 
     104           0 :   unsigned npoints = point_indices.size();
     105             :   // reserve space for the point vectors and values
     106           0 :   std::vector<std::vector<double>> points(npoints, std::vector<double>(dimension));
     107           0 :   std::vector<double> values(npoints);
     108             : 
     109             :   // fill point and value vectors for all the grid points
     110           0 :   for (unsigned i = 0; i < npoints; ++i) {
     111           0 :     points[i] = grid_pntr->getPoint(point_indices[i]);
     112           0 :     values[i] = grid_pntr->getValue(point_indices[i]);
     113             :   }
     114             : 
     115           0 :   return multiLinearInterpolation(arg, points, values, dimension);
     116           0 : }
     117             : 
     118             : 
     119     1263854 : double GridLinearInterpolation::getGridValueAndDerivativesWithLinearInterpolation_1D(GridBase* grid_pntr, const std::vector<double>& arg, std::vector<double>& der) {
     120     1263854 :   plumed_massert(grid_pntr->getDimension()==1,"The grid is of the wrong dimension, should be one-dimensional");
     121     1263854 :   plumed_massert(arg.size()==1,"input value is of the wrong size");
     122             : 
     123     1263854 :   double x = arg[0];
     124     1263854 :   double grid_dx = grid_pntr->getDx()[0];
     125     1263854 :   double grid_min; Tools::convert( grid_pntr->getMin()[0], grid_min);
     126             : 
     127     1263854 :   double xtoindex = (x-grid_min)/grid_dx;
     128     1263854 :   std::vector<unsigned int> i0(1); i0[0] = unsigned(std::floor(xtoindex));
     129     1263854 :   std::vector<unsigned int> i1(1); i1[0] = unsigned(std::ceil(xtoindex));
     130             :   //
     131     1263854 :   std::vector<double> d0 (1), d1 (1);
     132     1263854 :   double x0 = grid_pntr->getPoint(i0)[0];
     133     1263854 :   double x1 = grid_pntr->getPoint(i1)[0];
     134     1263854 :   double y0 = grid_pntr->getValueAndDerivatives(i0, d0);
     135     1263854 :   double y1 = grid_pntr->getValueAndDerivatives(i1, d1);
     136             :   //
     137     1263854 :   der.resize(1);
     138     2500429 :   der[0] = linearInterpolation(arg[0],x0,x1,d0[0],d1[0]);
     139     2527708 :   return linearInterpolation(arg[0],x0,x1,y0,y1);
     140             : }
     141             : 
     142             : 
     143           0 : double GridLinearInterpolation::getGridValueAndDerivativesWithLinearInterpolation_ND(GridBase* grid_pntr, const std::vector<double>& arg, std::vector<double>& der) {
     144           0 :   unsigned int dimension = grid_pntr->getDimension();
     145           0 :   plumed_massert(dimension==arg.size(),"The grid dimensions do not match the given arguments");
     146             :   // get points first
     147           0 :   std::vector<std::vector<unsigned>> point_indices = getAdjacentPoints(grid_pntr, arg);
     148             : 
     149           0 :   unsigned npoints = point_indices.size();
     150             :   // allocate space for the point vectors and values
     151           0 :   std::vector<std::vector<double>> points(npoints, std::vector<double>(dimension));
     152           0 :   std::vector<double> values(npoints);
     153           0 :   std::vector<std::vector<double>> derivs(npoints, std::vector<double>(dimension));
     154             : 
     155             :   // fill point, value and deriv vectors for all the grid points
     156           0 :   for (unsigned i = 0; i < npoints; ++i) {
     157           0 :     points[i] = grid_pntr->getPoint(point_indices[i]);
     158           0 :     values[i] = grid_pntr->getValueAndDerivatives(point_indices[i], derivs[i]);
     159             :   }
     160             : 
     161           0 :   for (size_t i = 0; i < derivs.size(); ++i) {
     162           0 :     der[i] = multiLinearInterpolation(arg, points, derivs[i], dimension);
     163             :   }
     164           0 :   return multiLinearInterpolation(arg, points, values, dimension);
     165           0 : }
     166             : 
     167             : 
     168       49033 : double GridLinearInterpolation::getGridValueWithLinearInterpolation(GridBase* grid_pntr, const std::vector<double>& arg) {
     169       49033 :   unsigned int dim = grid_pntr->getDimension();
     170       49033 :   if(dim==1) {
     171        1869 :     return getGridValueWithLinearInterpolation_1D(grid_pntr,arg);
     172             :   }
     173       47164 :   else if(dim==2) {
     174       47164 :     return getGridValueWithLinearInterpolation_2D(grid_pntr,arg);
     175             :   }
     176             :   else {
     177           0 :     return getGridValueWithLinearInterpolation_ND(grid_pntr,arg);
     178             :   }
     179             : }
     180             : 
     181             : 
     182     1263854 : double GridLinearInterpolation::getGridValueAndDerivativesWithLinearInterpolation(GridBase* grid_pntr, const std::vector<double>& arg, std::vector<double>& der) {
     183     1263854 :   unsigned int dim = grid_pntr->getDimension();
     184     1263854 :   if(dim==1) {
     185     1263854 :     return getGridValueAndDerivativesWithLinearInterpolation_1D(grid_pntr,arg,der);
     186             :   }
     187           0 :   return getGridValueAndDerivativesWithLinearInterpolation_ND(grid_pntr,arg,der);
     188             : }
     189             : 
     190             : 
     191             : // returns the adjacent Grid indices of all double arg as vector of vectors
     192           0 : std::vector<std::vector<unsigned>> GridLinearInterpolation::getAdjacentIndices(GridBase* grid_pntr, const std::vector<double>& arg) {
     193           0 :   unsigned int dimension = grid_pntr->getDimension();
     194             : 
     195           0 :   std::vector<std::vector<unsigned>> indices(dimension, std::vector<unsigned>(2));
     196           0 :   for (unsigned i=0; i<dimension; ++i) {
     197           0 :     std::vector<unsigned> temp_indices(2);
     198             :     //
     199           0 :     double grid_dx = grid_pntr->getDx()[i];
     200           0 :     double grid_min; Tools::convert( grid_pntr->getMin()[i], grid_min);
     201           0 :     double xtoindex = (arg[i]-grid_min)/grid_dx;
     202           0 :     temp_indices[0] = static_cast<unsigned>(std::floor(xtoindex));
     203           0 :     temp_indices[1] = static_cast<unsigned>(std::ceil(xtoindex));
     204           0 :     indices[i] = temp_indices;
     205             :   }
     206           0 :   return indices;
     207           0 : }
     208             : 
     209             : 
     210           0 : std::vector<std::vector<unsigned>> GridLinearInterpolation::getAdjacentPoints(GridBase* grid_pntr, const std::vector<double>& arg) {
     211             :   // upper and lower grid indices as vectors for each dimension
     212           0 :   std::vector<std::vector<unsigned>> grid_indices = getAdjacentIndices(grid_pntr, arg);
     213           0 :   unsigned npoints = 1U << grid_pntr->getDimension();
     214             : 
     215             :   // generate combination of grid indices if multidimensional to match the actual points
     216             :   // the retrieved combinations will be in column-major order
     217           0 :   std::vector<std::vector<unsigned>> point_indices(npoints);
     218           0 :   for (unsigned i = 0; i < npoints; ++i) {
     219             :     unsigned temp = i;
     220             :     std::vector<unsigned> current_indices;
     221           0 :     for (const auto& vec: grid_indices) {
     222           0 :       unsigned j = temp % 2;
     223           0 :       current_indices.push_back(vec[j]);
     224           0 :       temp /= 2;
     225             :     }
     226           0 :     point_indices[i] = current_indices;
     227             :   }
     228           0 :   return point_indices;
     229           0 : }
     230             : 
     231             : 
     232             : 
     233             : 
     234             : }
     235             : }