LCOV - code coverage report
Current view: top level - src - interpolation_engine.cxx (source / functions) Hit Total Coverage
Test: code_coverage_filter.info Lines: 1 257 0.4 %
Date: 2024-03-28 16:04:17 Functions: 2 16 12.5 %
Branches: 2 554 0.4 %

           Branch data     Line data    Source code
       1                 :            : /** 
       2                 :            :  *
       3                 :            :  *
       4                 :            :  * @author: William Kingsford (namespace created by Abed Yassine)
       5                 :            :  * @data: October 18th, 2019
       6                 :            :  *
       7                 :            :  * @file interpolation_engine.cxx
       8                 :            :  * @brief: contains implementations for the interpolation_engine methods
       9                 :            :  */
      10                 :            : 
      11                 :            : 
      12                 :            : #include "interpolation_engine.h"
      13                 :            : #include "profiler.h"
      14                 :            : 
      15                 :            : #include <cmath>
      16                 :            : 
      17                 :            : namespace SP_IE 
      18                 :            : {
      19                 :            : RaisedCosLUT *RaisedCosLUT::_s_instance = 0;
      20                 :            : 
      21                 :          0 : double sinc(double x) {
      22         [ #  # ]:          0 :     if (x == 0) {
      23                 :          0 :         return 1;
      24                 :            :     }
      25                 :            :     else {
      26                 :          0 :         return sin(M_PI * x) / (M_PI * x);
      27                 :            :     }
      28                 :            : }
      29                 :            : 
      30                 :          0 : double raised_cosine(double x, double beta) {
      31         [ #  # ]:          0 :     if (abs(x) == 1.0 / (2.0 * beta)) {
      32                 :          0 :         return (M_PI / (4.0)) * sinc(1/(2.0*beta));
      33                 :            :     }
      34                 :            :     else {
      35                 :          0 :         return sinc(x) * cos(M_PI*beta*x) / (1.0 - 4.0*beta*beta*x*x);
      36                 :            :     }
      37                 :            : }
      38                 :            : 
      39                 :          0 : double integrate_midpoint(std::function<double(double)> func, double l, double r, unsigned num_pts)
      40                 :            : {
      41                 :            :     /*
      42                 :            :      * This integration assumes that the delta_x is fixed. i.e, all points along the line are equidistant  
      43                 :            :      *
      44                 :            :      */
      45                 :            : 
      46                 :          0 :     double integration = 0.0;
      47                 :            :     
      48                 :          0 :     double temp_r = r; 
      49                 :          0 :     r = std::max(r, l); 
      50                 :          0 :     l = std::min(temp_r, l);
      51                 :            : 
      52         [ #  # ]:          0 :     if (num_pts == 0)
      53                 :            :     {
      54         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmintegrate_midpoint::num_points is zero ==> cannot integrate.\033[0;m\n", 31); 
      55                 :          0 :         std::exit(-1); 
      56                 :            :     }
      57                 :            : 
      58                 :          0 :     double delta_x = (r - l)/num_pts; 
      59                 :            :     
      60                 :          0 :     double x = l + 0.5*delta_x; 
      61         [ #  # ]:          0 :     for (unsigned i = 0; i < num_pts; i++)
      62                 :            :     {
      63         [ #  # ]:          0 :         integration += func(x); 
      64                 :          0 :         x += delta_x;
      65                 :            :     }
      66                 :            : 
      67                 :          0 :     return integration*delta_x;
      68                 :            : }
      69                 :            : 
      70                 :          0 : Eigen::RowVectorXd integrate_midpoint(std::function<Eigen::RowVectorXd(Eigen::RowVectorXd)> func,
      71                 :            :         Eigen::RowVectorXd l, Eigen::RowVectorXd r, unsigned num_pts)
      72                 :            : {
      73 [ #  # ][ #  # ]:          0 :     Eigen::RowVectorXd vec = r - l;
      74         [ #  # ]:          0 :     double delta_x = vec.norm()/num_pts; 
      75                 :            :     
      76 [ #  # ][ #  # ]:          0 :     std::cout << "l: " << l << std::endl;
                 [ #  # ]
      77         [ #  # ]:          0 :     vec.normalize();
      78                 :            :     
      79 [ #  # ][ #  # ]:          0 :     Eigen::RowVectorXd x = l + 0.5*delta_x*vec; 
                 [ #  # ]
      80                 :            :     
      81 [ #  # ][ #  # ]:          0 :     Eigen::RowVectorXd integration = func(x);
      82                 :            :     
      83         [ #  # ]:          0 :     for (unsigned i = 1; i < num_pts; i++)
      84                 :            :     {
      85 [ #  # ][ #  # ]:          0 :         std::cout << "x: " << x << std::endl;
                 [ #  # ]
      86 [ #  # ][ #  # ]:          0 :         x += delta_x*vec;
      87 [ #  # ][ #  # ]:          0 :         integration += func(x); 
                 [ #  # ]
      88                 :            :     }
      89                 :            : 
      90 [ #  # ][ #  # ]:          0 :     return integration*delta_x;
      91                 :            : }
      92                 :            : 
      93                 :          0 : double evaluate_interpolation(double x, double x_min, double x_max, double delta_x, 
      94                 :            :                               vector<double> &f, int N, bool use_raised_cosine, double beta) {
      95                 :            :     // determine nearest sample point index to the left of x
      96                 :          0 :     int idx = (int) floor((x - x_min)/delta_x);
      97                 :            :     // convolution will include N sample points on each side of x
      98                 :          0 :     int begin = idx - N + 1;
      99                 :          0 :     int end = idx + N;
     100                 :            : 
     101                 :            :     // ensure sum is within the range of the sample function
     102         [ #  # ]:          0 :     if (begin < 0) {
     103                 :          0 :         std::cout << "WARNING: Sum in evaluate_interpolation is being pruned " <<
     104                 :          0 :          " due to going outside of range (x = " << x << 
     105                 :          0 :          ", x_min = " << x_min << ", x_max = " << x_max << ")" << std::endl;
     106                 :          0 :         begin = 0;
     107                 :            :     }
     108         [ #  # ]:          0 :     if (end >= static_cast<int>(f.size())) {
     109                 :          0 :         std::cout << "WARNING: Sum in evaluate_interpolation is being " <<
     110                 :          0 :             " pruned due to going outside of range (x = " << x << 
     111                 :          0 :             ", x_min = " << x_min << ", x_max = " << x_max << ")" << std::endl;
     112                 :          0 :         end = static_cast<int>(f.size()) - 1;
     113                 :            :     }
     114                 :            :     
     115                 :          0 :     double ret_val = 0;
     116                 :          0 :     double a = (x - x_min)/delta_x; // pre-computing for efficiency
     117         [ #  # ]:          0 :     if (!use_raised_cosine) { // use sinc
     118         [ #  # ]:          0 :         for (int i=begin; i<=end; i++) {
     119                 :          0 :             ret_val += f[i] * sinc(a - (double)i);
     120                 :            :         }
     121                 :            :     }
     122                 :            :     else { // use raised cosine
     123         [ #  # ]:          0 :         for (int i=begin; i<=end; i++) {
     124                 :          0 :             ret_val += f[i] * raised_cosine(a - (double)i, beta);
     125                 :            :         }
     126                 :            :     }
     127                 :            : 
     128                 :          0 :     return ret_val;
     129                 :            : }
     130                 :            : 
     131                 :          0 : double evaluate_interpolation(double x[3], double x_min[3], double x_max[3], double delta_x, 
     132                 :            :                               vector<vector<vector<double>>> &f, int N, 
     133                 :            :                               bool use_raised_cosine, double beta) {
     134                 :            :     // determine nearest sample point index to the left of x
     135                 :            :     int idx[3], begin[3], end[3];
     136                 :            : 
     137                 :            :     // number of entries in f array across each dimension
     138                 :            :     int length[3];
     139                 :          0 :     length[0] = static_cast<int>(f.size());
     140                 :          0 :     length[1] = static_cast<int>(f[0].size());
     141                 :          0 :     length[2] = static_cast<int>(f[0][0].size());
     142         [ #  # ]:          0 :     for (int i=0; i<3; i++) {
     143                 :          0 :         idx[i] = (int) floor((x[i] - x_min[i])/delta_x);
     144                 :            : 
     145                 :            :         // convolution will include N sample points on each side of x
     146                 :          0 :         begin[i] = idx[i] - N + 1;
     147                 :          0 :         end[i] = idx[i] + N;
     148                 :            : 
     149                 :            :         // ensure sum is within the range of the sample function
     150         [ #  # ]:          0 :         if (begin[i] < 0) {
     151         [ #  # ]:          0 :             std::cout << "WARNING: Sum in evaluate_interpolation is being " << 
     152 [ #  # ][ #  # ]:          0 :                 "pruned due to going outside of range (x[" << i << "] = " << x[i] << 
         [ #  # ][ #  # ]
     153 [ #  # ][ #  # ]:          0 :                 ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     154                 :          0 :             begin[i] = 0;
     155                 :            :         }
     156         [ #  # ]:          0 :         if (end[i] >= static_cast<int>(f.size())) {
     157         [ #  # ]:          0 :             std::cout << "WARNING: Sum in evaluate_interpolation is being " <<
     158 [ #  # ][ #  # ]:          0 :                 "pruned due to going outside of range (x[" << i << "] = " << x[i] << 
         [ #  # ][ #  # ]
     159 [ #  # ][ #  # ]:          0 :                 ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     160                 :          0 :             end[i] = length[i] - 1;
     161                 :            :         }
     162                 :            :     }
     163                 :            :     
     164                 :          0 :     double ret_val = 0;
     165                 :            :     // pre-computing for efficiency
     166                 :          0 :     double a = (x[0] - x_min[0])/delta_x; 
     167                 :          0 :     double b = (x[1] - x_min[1])/delta_x; 
     168                 :          0 :     double c = (x[2] - x_min[2])/delta_x; 
     169                 :            : 
     170         [ #  # ]:          0 :     RaisedCosLUT *lut = RaisedCosLUT::instance();
     171                 :            :     // check if initialized already; if not, initialize with default settings
     172         [ #  # ]:          0 :     if (!lut->isInitialized()) {
     173         [ #  # ]:          0 :         lut->initialize((double)(N+1), 100.0, beta); // range of [-(N+1),N+1], resolution of 100
     174                 :            :     }
     175                 :            : 
     176         [ #  # ]:          0 :     if (!use_raised_cosine) { // use sinc
     177         [ #  # ]:          0 :         for (int i=begin[0]; i<=end[0]; i++) {
     178         [ #  # ]:          0 :             for (int j=begin[1]; j<=end[1]; j++) {
     179         [ #  # ]:          0 :                 for (int k=begin[2]; k<=end[2]; k++) {
     180         [ #  # ]:          0 :                     ret_val += f[i][j][k] * sinc(a - (double)i) * 
     181         [ #  # ]:          0 :                                             sinc(b - (double)j) *
     182         [ #  # ]:          0 :                                             sinc(c - (double)k);
     183                 :            :                 }
     184                 :            :             }
     185                 :            :         }
     186                 :            :     }
     187                 :            :     else { // use raised cosine
     188         [ #  # ]:          0 :         for (int i=begin[0]; i<=end[0]; i++) {
     189         [ #  # ]:          0 :             for (int j=begin[1]; j<=end[1]; j++) {
     190         [ #  # ]:          0 :                 for (int k=begin[2]; k<=end[2]; k++) {
     191                 :            :                     // we know that abs(a-(double)i) < N+1, so we don't need to worry about
     192                 :            :                     // going out of the array's bounds
     193                 :          0 :                     ret_val += f[i][j][k] * lut->get(a-(double)i) *
     194                 :          0 :                                             lut->get(b-(double)j) *
     195                 :          0 :                                             lut->get(c-(double)k);
     196                 :            :                     //ret_val += f[i][j][k] * raised_cosine(a-(double)i, beta) *
     197                 :            :                     //                        raised_cosine(b-(double)j, beta) *
     198                 :            :                     //                        raised_cosine(c-(double)k, beta);
     199                 :            :                 }
     200                 :            :             }
     201                 :            :         }
     202                 :            :     }
     203                 :            : 
     204                 :          0 :     return ret_val;
     205                 :            : }
     206                 :            : 
     207                 :            : // overloaded version for vector-valued functions of 2 variables
     208                 :          0 : Eigen::RowVectorXd evaluate_interpolation(double x[2], double x_min[2],
     209                 :            :                               double delta_x, 
     210                 :            :                               source_cache &grid_sources, int N, 
     211                 :            :                               bool use_raised_cosine, double beta) {
     212                 :            :     // determine nearest sample point index to the left of x
     213                 :            :     int idx[2], begin[2], end[2];
     214                 :            : 
     215                 :            :     
     216                 :            :     // number of entries in f array across each dimension
     217                 :            :     //int length[2];
     218                 :            :     //length[0] = f.size();
     219                 :            :     //length[1] = f[0].size();
     220         [ #  # ]:          0 :     for (int i=0; i<2; i++) {
     221                 :          0 :         idx[i] = (int) floor((x[i] - x_min[i])/delta_x);
     222                 :            : 
     223                 :            :         // convolution will include N sample points on each side of x
     224                 :          0 :         begin[i] = idx[i] - N + 1;
     225                 :          0 :         end[i] = idx[i] + N;
     226                 :            : 
     227                 :            :         // TODO: no longer need to check bounds now that we're using source_cache
     228                 :            :         // ensure sum is within the range of the sample function
     229                 :            :         /*if (begin[i] < 0) {
     230                 :            :             std::cout << "WARNING: Sum in evaluate_interpolation is being "
     231                 :            :                  "pruned due to going outside of range (x[" << i << "] = " << 
     232                 :            :                  x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
     233                 :            :             begin[i] = 0;
     234                 :            :         }
     235                 :            :         if (end[i] >= f.size()) {
     236                 :            :             std::cout << "WARNING: Sum in evaluate_interpolation is being "
     237                 :            :                 "pruned due to going outside of range (x[" << i << "] = " << 
     238                 :            :                 x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
     239                 :            :             end[i] = length[i] - 1;
     240                 :            :         }*/
     241                 :            :     }
     242                 :            :     
     243                 :            :     // initialize return value
     244 [ #  # ][ #  # ]:          0 :     Eigen::RowVectorXd ret_val = Eigen::RowVectorXd::Zero(grid_sources.get_num_tetras());
     245                 :            :     // pre-computing for efficiency
     246                 :          0 :     double a = (x[0] - x_min[0])/delta_x; 
     247                 :          0 :     double b = (x[1] - x_min[1])/delta_x; 
     248                 :            :     double weight;
     249                 :            : 
     250         [ #  # ]:          0 :     RaisedCosLUT *lut = RaisedCosLUT::instance();
     251                 :            :     // check if initialized already; if not, initialize with default settings
     252         [ #  # ]:          0 :     if (!lut->isInitialized()) {
     253                 :            :         //TODO: currently making larger than we'll ever need, should adjust this so it can
     254                 :            :         // check if the size is too small (for when N changes)
     255                 :            :         //lut->initialize((double)(N+1), 100.0, beta); // range of [-(N+1),N+1], resolution of 100
     256         [ #  # ]:          0 :         lut->initialize((double)(11), 100.0, beta); // range of [-(11),11], resolution of 100
     257                 :            :     }
     258                 :            : 
     259                 :            :     // TODO: using only linear interpolation for now
     260                 :            :     
     261                 :            :     // parameters for linear interpolation
     262                 :            :     double t_x, t_y; 
     263                 :            :     int n_x, n_y;
     264                 :          0 :     n_x = static_cast<int>(floor(a));
     265                 :          0 :     n_y = static_cast<int>(floor(b));
     266                 :          0 :     t_x = a - (double)n_x;
     267                 :          0 :     t_y = b - (double)n_y;
     268                 :            : 
     269                 :            :     //cout << "Evaluating linear interpolation: t_x = " << t_x << ", t_y = " << t_y << endl;
     270                 :            :     // TODO: no longer need to cast n_x, n_y to int
     271                 :            : 
     272         [ #  # ]:          0 :     if (!use_raised_cosine) { // use linear interpolation
     273         [ #  # ]:          0 :         ret_val = (1.0 - t_y) * (
     274 [ #  # ][ #  # ]:          0 :                        (1.0 - t_x) * grid_sources.get_source((int)n_x, (int)n_y) +
                 [ #  # ]
     275 [ #  # ][ #  # ]:          0 :                        t_x * grid_sources.get_source((int)n_x + 1, (int)n_y) ) +
                 [ #  # ]
     276         [ #  # ]:          0 :                   t_y * (
     277 [ #  # ][ #  # ]:          0 :                        (1.0 - t_x) * grid_sources.get_source((int)n_x, (int)n_y+1) +
                 [ #  # ]
     278 [ #  # ][ #  # ]:          0 :                        t_x * grid_sources.get_source((int)n_x + 1, (int)n_y+1) );
                 [ #  # ]
     279                 :            :         // TODO: old code for sinc
     280                 :            :         /*
     281                 :            :         for (int i=begin[0]; i<=end[0]; i++) {
     282                 :            :             for (int j=begin[1]; j<=end[1]; j++) {
     283                 :            :                 for (int k=begin[2]; k<=end[2]; k++) {
     284                 :            :                     weight = sinc(a - (double)i) * sinc(b - (double)j) * sinc(c - (double)k);
     285                 :            :                     ret_val += weight * f[i][j][k];
     286                 :            :                 }
     287                 :            :             }
     288                 :            :         }*/
     289                 :            :     }
     290                 :            :     else { // use raised cosine
     291         [ #  # ]:          0 :         for (int i=begin[0]; i<=end[0]; i++) {
     292         [ #  # ]:          0 :             for (int j=begin[1]; j<=end[1]; j++) {
     293                 :            :                 // we know that abs(a-(double)i) < N+1, so we don't need to worry about
     294                 :            :                 // going out of the array's bounds
     295                 :          0 :                 weight = lut->get(a-(double)i) *
     296                 :          0 :                          lut->get(b-(double)j);
     297                 :            :                 //ret_val += weight * f[i][j][k];
     298 [ #  # ][ #  # ]:          0 :                 ret_val += weight * grid_sources.get_source(i, j);
                 [ #  # ]
     299                 :            :             }
     300                 :            :         }
     301                 :            :     }
     302                 :            : 
     303                 :          0 :     return ret_val;
     304                 :            : }
     305                 :            : 
     306                 :            : // TODO: replace double[2] with std::vector
     307                 :          0 : vector<Eigen::RowVectorXd> evaluate_interpolation_gradient(double x[2], double x_min[2],
     308                 :            :                               double delta_x, 
     309                 :            :                               source_cache &grid_sources, int N, 
     310                 :            :                               bool /*use_raised_cosine*/, double beta) {
     311                 :            :     // TODO: don't need the below now that we're using source_cache
     312                 :            :     /*
     313                 :            :     // determine nearest sample point index to the left of x
     314                 :            :     int idx[2], begin[2], end[2];
     315                 :            : 
     316                 :            :     // number of entries in f array across each dimension
     317                 :            :     int length[2];
     318                 :            :     length[0] = f.size();
     319                 :            :     length[1] = f[0].size();
     320                 :            :     for (int i=0; i<2; i++) {
     321                 :            :         idx[i] = (int) floor((x[i] - x_min[i])/delta_x);
     322                 :            : 
     323                 :            :         // convolution will include N sample points on each side of x
     324                 :            :         begin[i] = idx[i] - N + 1;
     325                 :            :         end[i] = idx[i] + N;
     326                 :            : 
     327                 :            :         // ensure sum is within the range of the sample function
     328                 :            :         if (begin[i] < 0) {
     329                 :            :             std::cout << "WARNING: Sum in evaluate_interpolation is being " << 
     330                 :            :                 "pruned due to going outside of range (x[" << i << "] = " << 
     331                 :            :                 x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
     332                 :            :             begin[i] = 0;
     333                 :            :         }
     334                 :            :         if (end[i] >= f.size()) {
     335                 :            :             std::cout << "WARNING: Sum in evaluate_interpolation is being " << 
     336                 :            :                 "pruned due to going outside of range (x[" << i << "] = " << 
     337                 :            :                 x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
     338                 :            :             end[i] = length[i] - 1;
     339                 :            :         }
     340                 :            :     }*/
     341                 :            :     
     342                 :            :     // initialize return value
     343                 :            :     // TODO: replace 2 by the number of components in f
     344 [ #  # ][ #  # ]:          0 :     vector<Eigen::RowVectorXd> ret_val(2, Eigen::RowVectorXd::Zero(grid_sources.get_num_tetras()));
                 [ #  # ]
     345                 :            :     // pre-computing for efficiency
     346                 :          0 :     double a = (x[0] - x_min[0])/delta_x; 
     347                 :          0 :     double b = (x[1] - x_min[1])/delta_x; 
     348                 :            : 
     349         [ #  # ]:          0 :     RaisedCosLUT *lut = RaisedCosLUT::instance();
     350                 :            :     // check if initialized already; if not, initialize with default settings
     351         [ #  # ]:          0 :     if (!lut->isInitialized()) {
     352         [ #  # ]:          0 :         lut->initialize((double)(N+1), 100.0, beta); // range of [-(N+1),N+1], resolution of 100
     353                 :            :     }
     354                 :            : 
     355                 :            :     // TODO: using only linear interpolation for now
     356                 :            :     
     357                 :            :     // parameters for linear interpolation
     358                 :            :     double t_x, t_y; 
     359                 :            :     int n_x, n_y;
     360                 :          0 :     n_x = static_cast<int>(floor(a));
     361                 :          0 :     n_y = static_cast<int>(floor(b));
     362                 :          0 :     t_x = a - (double)n_x;
     363                 :          0 :     t_y = b - (double)n_y;
     364                 :            : 
     365                 :            :     //cout << "Evaluating linear interpolation: t_x = " << t_x << ", t_y = " << t_y << endl;
     366                 :            :     // TODO: use intermediate variables to make this less ugly (e.g. variables so we don't call
     367                 :            :     // grid_sources.get_source so often)
     368 [ #  # ][ #  # ]:          0 :     ret_val[0] = (1.0 / delta_x) * ( (1.0 - t_y) * (
     369 [ #  # ][ #  # ]:          0 :                 grid_sources.get_source(n_x + 1, n_y) - grid_sources.get_source(n_x, n_y) ) +
         [ #  # ][ #  # ]
     370         [ #  # ]:          0 :                  t_y * (
     371 [ #  # ][ #  # ]:          0 :                 grid_sources.get_source(n_x + 1, n_y+1) - grid_sources.get_source(n_x, n_y+1) )
                 [ #  # ]
     372         [ #  # ]:          0 :                            );
     373         [ #  # ]:          0 :     ret_val[1] = (1.0 / delta_x) * ( 
     374 [ #  # ][ #  # ]:          0 :                   ( (1.0 - t_x) * grid_sources.get_source(n_x, n_y+1) +
                 [ #  # ]
     375 [ #  # ][ #  # ]:          0 :                    t_x * grid_sources.get_source(n_x + 1, n_y+1) )
     376         [ #  # ]:          0 :                   -
     377 [ #  # ][ #  # ]:          0 :                   ( (1.0 - t_x) * grid_sources.get_source(n_x, n_y) +
                 [ #  # ]
     378 [ #  # ][ #  # ]:          0 :                    t_x * grid_sources.get_source(n_x + 1, n_y) )
     379         [ #  # ]:          0 :                      );
     380                 :            :     /*ret_val[0] = (1.0 / delta_x) * ( (1.0 - t_y) * (
     381                 :            :                                f[n_x + 1][n_y] - f[n_x][n_y] ) +
     382                 :            :                              t_y * (
     383                 :            :                                f[n_x + 1][n_y+1] - f[n_x][n_y+1] )
     384                 :            :                            );
     385                 :            :     ret_val[1] = (1.0 / delta_x) * ( 
     386                 :            :                   ( (1.0 - t_x) * f[n_x][n_y+1] +
     387                 :            :                    t_x * f[n_x + 1][n_y+1] )
     388                 :            :                   -
     389                 :            :                   ( (1.0 - t_x) * f[n_x][n_y] +
     390                 :            :                    t_x * f[n_x + 1][n_y] )
     391                 :            :                      );*/
     392                 :            : 
     393                 :            :     // TODO: uncomment, add an option for linear interpolation
     394                 :            :     /*
     395                 :            :     double weight;
     396                 :            :     if (!use_raised_cosine) { // use sinc
     397                 :            :         for (int i=begin[0]; i<=end[0]; i++) {
     398                 :            :             for (int j=begin[1]; j<=end[1]; j++) {
     399                 :            :                 for (int k=begin[2]; k<=end[2]; k++) {
     400                 :            :                     weight = sinc(a - (double)i) * sinc(b - (double)j) * sinc(c - (double)k);
     401                 :            :                     ret_val += weight * f[i][j][k];
     402                 :            :                 }
     403                 :            :             }
     404                 :            :         }
     405                 :            :     }
     406                 :            :     else { // use raised cosine
     407                 :            :         for (int i=begin[0]; i<=end[0]; i++) {
     408                 :            :             for (int j=begin[1]; j<=end[1]; j++) {
     409                 :            :                 for (int k=begin[2]; k<=end[2]; k++) {
     410                 :            :                     // we know that abs(a-(double)i) < N+1, so we don't need to worry about
     411                 :            :                     // going out of the array's bounds
     412                 :            :                     weight = lut->get(a-(double)i) *
     413                 :            :                              lut->get(b-(double)j) *
     414                 :            :                              lut->get(c-(double)k);
     415                 :            :                     ret_val += weight * f[i][j][k];
     416                 :            :                 }
     417                 :            :             }
     418                 :            :         }
     419                 :            :     }*/
     420                 :            : 
     421                 :          0 :     return ret_val;
     422                 :            : }
     423                 :            : 
     424                 :            : // overloaded version for vector-valued functions of 3 variables
     425                 :          0 : Eigen::RowVectorXd evaluate_interpolation(double x[3], double x_min[3], double x_max[3],
     426                 :            :                               double delta_x, 
     427                 :            :                               vector<vector<vector<Eigen::RowVectorXd>>> &f, int N, 
     428                 :            :                               bool /*use_raised_cosine*/, double beta) {
     429                 :            :     // determine nearest sample point index to the left of x
     430                 :            :     int idx[3], begin[3], end[3];
     431                 :            : 
     432                 :            :     // number of entries in f array across each dimension
     433                 :            :     int length[3];
     434                 :          0 :     length[0] = static_cast<int>(f.size());
     435                 :          0 :     length[1] = static_cast<int>(f[0].size());
     436                 :          0 :     length[2] = static_cast<int>(f[0][0].size());
     437         [ #  # ]:          0 :     for (int i=0; i<3; i++) {
     438                 :          0 :         idx[i] = (int) floor((x[i] - x_min[i])/delta_x);
     439                 :            : 
     440                 :            :         // convolution will include N sample points on each side of x
     441                 :          0 :         begin[i] = idx[i] - N + 1;
     442                 :          0 :         end[i] = idx[i] + N;
     443                 :            : 
     444                 :            :         // ensure sum is within the range of the sample function
     445         [ #  # ]:          0 :         if (begin[i] < 0) {
     446         [ #  # ]:          0 :             std::cout << "WARNING: Sum in evaluate_interpolation is being " << 
     447 [ #  # ][ #  # ]:          0 :                 "pruned due to going outside of range (x[" << i << "] = " << 
                 [ #  # ]
     448 [ #  # ][ #  # ]:          0 :                 x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     449                 :          0 :             begin[i] = 0;
     450                 :            :         }
     451         [ #  # ]:          0 :         if (end[i] >= static_cast<int>(f.size())) {
     452         [ #  # ]:          0 :             std::cout << "WARNING: Sum in evaluate_interpolation is being " << 
     453 [ #  # ][ #  # ]:          0 :                 "pruned due to going outside of range (x[" << i << "] = " << 
                 [ #  # ]
     454 [ #  # ][ #  # ]:          0 :                 x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     455                 :          0 :             end[i] = length[i] - 1;
     456                 :            :         }
     457                 :            :     }
     458                 :            :     
     459                 :            :     // initialize return value
     460 [ #  # ][ #  # ]:          0 :     Eigen::RowVectorXd ret_val = Eigen::RowVectorXd::Zero(f[0][0][0].size());
                 [ #  # ]
     461                 :            :     // pre-computing for efficiency
     462                 :          0 :     double a = (x[0] - x_min[0])/delta_x; 
     463                 :          0 :     double b = (x[1] - x_min[1])/delta_x; 
     464                 :          0 :     double c = (x[2] - x_min[2])/delta_x; 
     465                 :            : 
     466         [ #  # ]:          0 :     RaisedCosLUT *lut = RaisedCosLUT::instance();
     467                 :            :     // check if initialized already; if not, initialize with default settings
     468         [ #  # ]:          0 :     if (!lut->isInitialized()) {
     469         [ #  # ]:          0 :         lut->initialize((double)(N+1), 100.0, beta); // range of [-(N+1),N+1], resolution of 100
     470                 :            :     }
     471                 :            : 
     472                 :            :     // TODO: using only linear interpolation for now
     473                 :            :     
     474                 :            :     // parameters for linear interpolation
     475                 :            :     double t_x, t_y, t_z; 
     476                 :            :     double n_x, n_y, n_z;
     477                 :          0 :     n_x = floor(a);
     478                 :          0 :     n_y = floor(b);
     479                 :          0 :     n_z = floor(c);
     480                 :          0 :     t_x = a - n_x;
     481                 :          0 :     t_y = b - n_y;
     482                 :          0 :     t_z = c - n_z;
     483                 :            : 
     484 [ #  # ][ #  # ]:          0 :     cout << "Evaluating linear interpolation..." << endl;
     485         [ #  # ]:          0 :     ret_val = (1.0 - t_z) * (
     486         [ #  # ]:          0 :                   (1.0 - t_y) * (
     487 [ #  # ][ #  # ]:          0 :                        (1.0 - t_x) * f[(int)n_x][(int)n_y][(int)n_z] +
     488 [ #  # ][ #  # ]:          0 :                        t_x * f[(int)n_x + 1][(int)n_y][(int)n_z] ) +
     489         [ #  # ]:          0 :                   t_y * (
     490 [ #  # ][ #  # ]:          0 :                        (1.0 - t_x) * f[(int)n_x][(int)n_y+1][(int)n_z] +
     491 [ #  # ][ #  # ]:          0 :                        t_x * f[(int)n_x + 1][(int)n_y+1][(int)n_z] ) ) +
     492         [ #  # ]:          0 :               t_z * (
     493         [ #  # ]:          0 :                   (1.0 - t_y) * (
     494 [ #  # ][ #  # ]:          0 :                        (1.0 - t_x) * f[(int)n_x][(int)n_y][(int)n_z+1] +
     495 [ #  # ][ #  # ]:          0 :                        t_x * f[(int)n_x + 1][(int)n_y][(int)n_z+1] ) +
     496         [ #  # ]:          0 :                   t_y * (
     497 [ #  # ][ #  # ]:          0 :                        (1.0 - t_x) * f[(int)n_x][(int)n_y+1][(int)n_z+1] +
     498 [ #  # ][ #  # ]:          0 :                        t_x * f[(int)n_x + 1][(int)n_y+1][(int)n_z+1] ) );
     499                 :            : 
     500                 :            :     // TODO: uncomment, add an option for linear interpolation
     501                 :            :     /*
     502                 :            :     double weight;
     503                 :            :     if (!use_raised_cosine) { // use sinc
     504                 :            :         for (int i=begin[0]; i<=end[0]; i++) {
     505                 :            :             for (int j=begin[1]; j<=end[1]; j++) {
     506                 :            :                 for (int k=begin[2]; k<=end[2]; k++) {
     507                 :            :                     weight = sinc(a - (double)i) * sinc(b - (double)j) * sinc(c - (double)k);
     508                 :            :                     ret_val += weight * f[i][j][k];
     509                 :            :                 }
     510                 :            :             }
     511                 :            :         }
     512                 :            :     }
     513                 :            :     else { // use raised cosine
     514                 :            :         for (int i=begin[0]; i<=end[0]; i++) {
     515                 :            :             for (int j=begin[1]; j<=end[1]; j++) {
     516                 :            :                 for (int k=begin[2]; k<=end[2]; k++) {
     517                 :            :                     // we know that abs(a-(double)i) < N+1, so we don't need to worry about
     518                 :            :                     // going out of the array's bounds
     519                 :            :                     weight = lut->get(a-(double)i) *
     520                 :            :                              lut->get(b-(double)j) *
     521                 :            :                              lut->get(c-(double)k);
     522                 :            :                     ret_val += weight * f[i][j][k];
     523                 :            :                 }
     524                 :            :             }
     525                 :            :         }
     526                 :            :     }*/
     527                 :            : 
     528                 :          0 :     return ret_val;
     529                 :            : }
     530                 :            : 
     531                 :          0 : Eigen::RowVectorXd evaluate_interpolation(const Eigen::RowVectorXd & pt, 
     532                 :            :                                 const Eigen::RowVectorXd & x_min, 
     533                 :            :                                 double delta_x,
     534                 :            :                                 source_cache& grid_sources)
     535                 :            : {
     536                 :            :     // initialize return value
     537 [ #  # ][ #  # ]:          0 :     Eigen::RowVectorXd ret_val = Eigen::RowVectorXd::Zero(grid_sources.get_num_tetras());
     538                 :            :     // pre-computing for efficiency
     539 [ #  # ][ #  # ]:          0 :     double a = (pt[0] - x_min[0])/delta_x; 
     540 [ #  # ][ #  # ]:          0 :     double b = (pt[1] - x_min[1])/delta_x; 
     541 [ #  # ][ #  # ]:          0 :     double c = (pt[2] - x_min[2])/delta_x; 
     542                 :            : 
     543                 :            :     // parameters for linear interpolation
     544                 :            :     double t_x, t_y, t_z; 
     545                 :            :     int n_x, n_y, n_z;
     546                 :          0 :     n_x = static_cast<int>(floor(a));
     547                 :          0 :     n_y = static_cast<int>(floor(b));
     548                 :          0 :     n_z = static_cast<int>(floor(c));
     549                 :          0 :     t_x = a - (double)n_x;
     550                 :          0 :     t_y = b - (double)n_y;
     551                 :          0 :     t_z = c - (double)n_z;
     552                 :            : 
     553 [ #  # ][ #  # ]:          0 :     cout << "Evaluating linear interpolation..." << endl;
     554         [ #  # ]:          0 :     ret_val = (1.0 - t_z) * (
     555         [ #  # ]:          0 :                   (1.0 - t_y) * (
     556 [ #  # ][ #  # ]:          0 :                        (1.0 - t_x) * grid_sources.get_source(n_x, n_y, n_z) +
                 [ #  # ]
     557 [ #  # ][ #  # ]:          0 :                        t_x * grid_sources.get_source(n_x + 1, n_y, n_z) ) +
                 [ #  # ]
     558         [ #  # ]:          0 :                   t_y * (
     559 [ #  # ][ #  # ]:          0 :                        (1.0 - t_x) * grid_sources.get_source(n_x, n_y+1, n_z) +
                 [ #  # ]
     560 [ #  # ][ #  # ]:          0 :                        t_x * grid_sources.get_source(n_x + 1, n_y+1, n_z) ) ) +
                 [ #  # ]
     561         [ #  # ]:          0 :               t_z * (
     562         [ #  # ]:          0 :                   (1.0 - t_y) * (
     563 [ #  # ][ #  # ]:          0 :                        (1.0 - t_x) * grid_sources.get_source(n_x, n_y, n_z+1) +
                 [ #  # ]
     564 [ #  # ][ #  # ]:          0 :                        t_x * grid_sources.get_source(n_x + 1, n_y, n_z+1) ) +
                 [ #  # ]
     565         [ #  # ]:          0 :                   t_y * (
     566 [ #  # ][ #  # ]:          0 :                        (1.0 - t_x) * grid_sources.get_source(n_x, n_y+1, n_z+1) +
                 [ #  # ]
     567 [ #  # ][ #  # ]:          0 :                        t_x * grid_sources.get_source(n_x + 1, n_y+1, n_z+1) ) );
                 [ #  # ]
     568                 :            : 
     569                 :          0 :     return ret_val;
     570                 :            : }
     571                 :            : 
     572                 :          0 : void sample_function(std::function<double(double)> &f,
     573                 :            :                                double x_min,
     574                 :            :                                double x_max,
     575                 :            :                                double step,
     576                 :            :                                vector<double> & func_values) {
     577                 :            :     // allocate memory
     578                 :          0 :     int n_samples = (int)floor((x_max - x_min)/step);
     579                 :          0 :     func_values.resize(n_samples);
     580                 :            : 
     581                 :            :     // sample function
     582         [ #  # ]:          0 :     for (int i = 0; i<n_samples; i++) {
     583                 :          0 :         func_values[i] = f(x_min + step*i);
     584                 :            :     }
     585                 :          0 : }
     586                 :            : 
     587                 :            : // overloaded 3d version
     588                 :          0 : void sample_function(std::function<double(double, double, double)> &f,
     589                 :            :                                double x_min[3],
     590                 :            :                                double x_max[3],
     591                 :            :                                double step,
     592                 :            :                                vector<vector<vector<double>>> & func_values) {
     593                 :            :     // allocate memory
     594                 :            :     int n_samples[3];
     595         [ #  # ]:          0 :     for (int i=0; i<3; i++) {
     596                 :          0 :         n_samples[i] = (int)floor((x_max[i] - x_min[i])/step);
     597                 :            :     }
     598 [ #  # ][ #  # ]:          0 :     func_values.resize(n_samples[0], vector<vector<double>>(n_samples[1],
     599         [ #  # ]:          0 :                                                             vector<double>(n_samples[2])));
     600                 :            :     // sample function
     601         [ #  # ]:          0 :     for (int i = 0; i<n_samples[0]; i++) {
     602         [ #  # ]:          0 :         for (int j = 0; j<n_samples[1]; j++) {
     603         [ #  # ]:          0 :             for (int k = 0; k<n_samples[2]; k++) {
     604         [ #  # ]:          0 :                 func_values[i][j][k] = f(x_min[0] + step*i, x_min[1] + step*j, x_min[2] + step*k);
     605                 :            :             }
     606                 :            :         }
     607                 :            :     }
     608                 :          0 : }
     609                 :            : 
     610                 :          0 : vector<float> interpolate_line_source(const SP_Line& line, const vector<double> & x_min, 
     611                 :            :                     source_cache & grid_sources, unsigned num_pts, double delta_x)
     612                 :            : {
     613                 :            :     // first define the function that would compute the fluence from a point source
     614                 :          0 :     std::function<Eigen::RowVectorXd(Eigen::RowVectorXd)> fluence_pt = [&](Eigen::RowVectorXd t)->Eigen::RowVectorXd 
     615                 :            :     {
     616                 :            :         // 1- define delta_x TODO
     617         [ #  # ]:          0 :         Eigen::RowVectorXd xmin(3); 
     618         [ #  # ]:          0 :         for (unsigned i = 0; i < 3; i++) {
     619         [ #  # ]:          0 :             xmin(i) = x_min[i]; 
     620                 :            :         }
     621                 :            :         // 2- Call evaluate interpolation for point sources that make use of grid_sources
     622         [ #  # ]:          0 :         return SP_IE::evaluate_interpolation(t, xmin, delta_x, grid_sources); 
     623         [ #  # ]:          0 :     };
     624                 :            : 
     625 [ #  # ][ #  # ]:          0 :     Eigen::RowVectorXd l(3), r(3); 
     626 [ #  # ][ #  # ]:          0 :     l(0) = static_cast<double>(line.end0.get_xcoord());
     627 [ #  # ][ #  # ]:          0 :     l(1) = static_cast<double>(line.end0.get_ycoord());
     628 [ #  # ][ #  # ]:          0 :     l(2) = static_cast<double>(line.end0.get_zcoord());
     629 [ #  # ][ #  # ]:          0 :     r(0) = static_cast<double>(line.end1.get_xcoord());
     630 [ #  # ][ #  # ]:          0 :     r(1) = static_cast<double>(line.end1.get_ycoord());
     631 [ #  # ][ #  # ]:          0 :     r(2) = static_cast<double>(line.end1.get_zcoord());
     632 [ #  # ][ #  # ]:          0 :     Eigen::RowVectorXd fluence_vals = integrate_midpoint(fluence_pt, l, r, num_pts); 
         [ #  # ][ #  # ]
     633                 :            : 
     634                 :            :     // TODO: might be expensive
     635                 :          0 :     vector<float> fluence;
     636 [ #  # ][ #  # ]:          0 :     for (unsigned i = 0; i < fluence_vals.size(); i++) {
     637 [ #  # ][ #  # ]:          0 :         fluence.push_back(static_cast<float>(fluence_vals(i))); 
     638                 :            :     }
     639                 :            : 
     640                 :          0 :     return fluence;
     641                 :            : }
     642                 :            : 
     643 [ +  - ][ +  - ]:         44 : };

Generated by: LCOV version 1.12