Line data    Source code

       1             : /****************************************************************************/
       2             : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
       3             : // Copyright (C) 2002-2024 German Aerospace Center (DLR) and others.
       4             : // This program and the accompanying materials are made available under the
       5             : // terms of the Eclipse Public License 2.0 which is available at
       6             : // https://www.eclipse.org/legal/epl-2.0/
       7             : // This Source Code may also be made available under the following Secondary
       8             : // Licenses when the conditions for such availability set forth in the Eclipse
       9             : // Public License 2.0 are satisfied: GNU General Public License, version 2
      10             : // or later which is available at
      11             : // https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
      12             : // SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
      13             : /****************************************************************************/
      14             : /// @file    CharacteristicMap.cpp
      15             : /// @author  Kevin Badalian (badalian_k@mmp.rwth-aachen.de)
      16             : /// @date    2021-02
      17             : ///
      18             : // Characteristic map for vehicle type parameters as needed by the MMPEVEM model
      19             : // Teaching and Research Area Mechatronics in Mobile Propulsion (MMP), RWTH Aachen
      20             : /****************************************************************************/
      21             : #include <config.h>
      22             : 
      23             : #include <cmath>
      24             : #include <cstring>
      25             : #include <stdexcept>
      26             : 
      27             : #include <utils/common/StringTokenizer.h>
      28             : #include <utils/emissions/CharacteristicMap.h>
      29             : 
      30             : 
      31             : // ===========================================================================
      32             : // method definitions
      33             : // ===========================================================================
      34             : void
      35      276162 : CharacteristicMap::determineStrides() {
      36      276162 :     strides.clear();
      37      276162 :     strides.reserve(domainDim);
      38      276162 :     strides.push_back(imageDim);
      39      552324 :     for (int i = 1; i < domainDim; i++) {
      40      276162 :         strides.push_back((int)axes[i - 1].size() * strides[i - 1]);
      41             :     }
      42      276162 : }
      43             : 
      44             : 
      45             : int
      46        9740 : CharacteristicMap::calcFlatIdx(const std::vector<int>& ref_idxs) const {
      47        9740 :     if (static_cast<int>(ref_idxs.size()) != domainDim) {
      48           0 :         throw std::runtime_error("The number of indices differs from the map's"
      49           0 :                                  " domain dimension.");
      50             :     }
      51             : 
      52             :     int flatIdx = 0;
      53       29220 :     for (int i = 0; i < domainDim; i++) {
      54       19480 :         if (ref_idxs[i] < 0) {
      55           0 :             throw std::runtime_error("The argument indices aren't non-negative.");
      56             :         }
      57       19480 :         flatIdx += ref_idxs[i] * strides[i];
      58             :     }
      59        9740 :     return flatIdx;
      60             : }
      61             : 
      62             : 
      63             : int
      64        5280 : CharacteristicMap::findNearestNeighborIdxs(const std::vector<double>& ref_p,
      65             :         std::vector<int>& ref_idxs, double eps) const {
      66        5280 :     if (static_cast<int>(ref_p.size()) != domainDim) {
      67           0 :         throw std::runtime_error("The argument point's size doesn't match the"
      68           0 :                                  " domain dimension.");
      69             :     }
      70             : 
      71        5280 :     ref_idxs = std::vector<int>(domainDim, -1);
      72       15840 :     for (int i = 0; i < domainDim; i++) {
      73       10560 :         if (axes[i][0] - eps <= ref_p[i] && ref_p[i] < axes[i][0]) {
      74           0 :             ref_idxs[i] = 0;
      75       10560 :         } else if (axes[i][axes[i].size() - 1] <= ref_p[i]
      76       10560 :                    && ref_p[i] < axes[i][axes[i].size() - 1] + eps) {
      77           0 :             ref_idxs[i] = (int)axes[i].size() - 1;
      78             :         } else {
      79       73418 :             for (int j = 0; j < static_cast<int>(axes[i].size()) - 1; j++) {
      80       73418 :                 if (axes[i][j] <= ref_p[i] && ref_p[i] < axes[i][j + 1]) {
      81             :                     // Pick the index that is closest to the point
      82       10560 :                     if (ref_p[i] - axes[i][j] <= axes[i][j + 1] - ref_p[i]) {
      83        8252 :                         ref_idxs[i] = j;
      84        8252 :                         break;
      85             :                     } else {
      86        2308 :                         ref_idxs[i] = j + 1;
      87        2308 :                         break;
      88             :                     }
      89             :                 }
      90             :             }
      91             :         }
      92             : 
      93             :         // The point lies outside of the valid range
      94       10560 :         if (ref_idxs[i] == -1) {
      95             :             return -1;
      96             :         }
      97             :     }
      98             : 
      99             :     return 0;
     100             : }
     101             : 
     102             : 
     103             : std::vector<double>
     104        9740 : CharacteristicMap::at(const std::vector<int>& ref_idxs) const {
     105        9740 :     if (static_cast<int>(ref_idxs.size()) != domainDim) {
     106           0 :         throw std::runtime_error("The number of indices differs from the map's"
     107           0 :                                  " domain dimension.");
     108             :     }
     109             : 
     110        9740 :     int flatIdx = calcFlatIdx(ref_idxs);
     111             :     return std::vector<double>(flattenedMap.begin() + flatIdx,
     112        9740 :                                flattenedMap.begin() + flatIdx + imageDim);
     113             : }
     114             : 
     115             : 
     116           0 : CharacteristicMap::CharacteristicMap(int domainDim, int imageDim,
     117             :                                      const std::vector<std::vector<double>>& ref_axes,
     118           0 :                                      const std::vector<double>& ref_flattenedMap)
     119           0 :     : domainDim(domainDim),
     120           0 :       imageDim(imageDim),
     121           0 :       axes(ref_axes),
     122           0 :       flattenedMap(ref_flattenedMap) {
     123             :     // Check whether the dimensions are consistent
     124           0 :     if (static_cast<int>(axes.size()) != domainDim) {
     125           0 :         throw std::runtime_error("The number of axes doesn't match the specified"
     126           0 :                                  " domain dimension.");
     127             :     }
     128             :     int expectedEntryCnt = imageDim;
     129           0 :     for (auto& ref_axis : axes) {
     130           0 :         expectedEntryCnt *= (int)ref_axis.size();
     131             :     }
     132           0 :     if (static_cast<int>(flattenedMap.size()) != expectedEntryCnt) {
     133           0 :         throw std::runtime_error("The number of map entries isn't equal to the"
     134           0 :                                  " product of the axes' dimensions times the image dimension.");
     135             :     }
     136             : 
     137           0 :     determineStrides();
     138           0 : }
     139             : 
     140             : 
     141      276162 : CharacteristicMap::CharacteristicMap(const std::string& ref_mapString) {
     142             :     // Split the map string into its three main parts
     143     1104648 :     const std::vector<std::string> tokens = StringTokenizer(ref_mapString, "|").getVector();
     144      276162 :     if (tokens.size() != 3) {
     145           0 :         throw std::runtime_error("The map string isn't made up of the 3 parts"
     146           0 :                                  " dimensions, axes, and flattened entries.");
     147             :     }
     148             : 
     149             :     // Extract the domain and image dimensions
     150      828486 :     const std::vector<std::string> dimensionTokens = StringTokenizer(tokens[0], ",").getVector();
     151      276162 :     if (dimensionTokens.size() != 2) {
     152           0 :         throw std::runtime_error("The domain and image dimensions aren't specified"
     153           0 :                                  " correctly.");
     154             :     }
     155      276162 :     domainDim = std::stoi(dimensionTokens[0]);
     156      276162 :     imageDim = std::stoi(dimensionTokens[1]);
     157             : 
     158             :     // Create the map axes
     159     1104648 :     const std::vector<std::string> axisTokens = StringTokenizer(tokens[1], ";").getVector();
     160      276162 :     if (static_cast<int>(axisTokens.size()) != domainDim) {
     161           0 :         throw std::runtime_error("The number of axes doesn't match the specified"
     162           0 :                                  " domain dimension.");
     163             :     }
     164      828486 :     for (auto& ref_axisToken : axisTokens) {
     165     1657032 :         std::vector<std::string> axisEntryTokens = StringTokenizer(ref_axisToken, ",").getVector();
     166             :         std::vector<double> axisEntries;
     167     1658748 :         for (auto& ref_axisEntryToken : axisEntryTokens) {
     168     1106424 :             axisEntries.push_back(std::stod(ref_axisEntryToken));
     169             :         }
     170      552324 :         axes.push_back(axisEntries);
     171      552324 :     }
     172             : 
     173             :     // Create the flattened map
     174      828516 :     const std::vector<std::string> flattenedMapTokens = StringTokenizer(tokens[2], ",").getVector();
     175      276162 :     int expectedEntryCnt = imageDim;
     176      828486 :     for (auto& ref_axis : axes) {
     177      552324 :         expectedEntryCnt *= (int)ref_axis.size();
     178             :     }
     179      276162 :     if (static_cast<int>(flattenedMapTokens.size()) != expectedEntryCnt) {
     180           0 :         throw std::runtime_error("The number of map entries isn't equal to the"
     181           0 :                                  " product of the axes' dimensions times the image dimension.");
     182             :     }
     183      276162 :     flattenedMap.reserve(expectedEntryCnt);
     184     1410570 :     for (auto& ref_flattenedMapToken : flattenedMapTokens) {
     185     1134408 :         flattenedMap.push_back(std::stod(ref_flattenedMapToken));
     186             :     }
     187             : 
     188      276162 :     determineStrides();
     189      276162 : }
     190             : 
     191             : 
     192             : std::string
     193           0 : CharacteristicMap::toString() const {
     194             :     // Write the domain and image dimensions
     195           0 :     std::string mapString = std::to_string(domainDim) + ","
     196           0 :                             + std::to_string(imageDim) + "|";
     197             : 
     198             :     // Add the axes
     199           0 :     for (int i = 0; i < static_cast<int>(axes.size()); i++) {
     200           0 :         for (int j = 0; j < static_cast<int>(axes[i].size()); j++) {
     201           0 :             mapString += std::to_string(axes[i][j])
     202           0 :                          + (j == static_cast<int>(axes[i].size()) - 1 ? "" : ",");
     203             :         }
     204           0 :         mapString += (i == static_cast<int>(axes.size()) - 1 ? "|" : ";");
     205             :     }
     206             : 
     207             :     // Append the flattened map entries
     208           0 :     for (int i = 0; i < static_cast<int>(flattenedMap.size()); i++) {
     209           0 :         mapString += std::to_string(flattenedMap[i])
     210           0 :                      + (i == static_cast<int>(flattenedMap.size()) - 1 ? "" : ",");
     211             :     }
     212             : 
     213           0 :     return mapString;
     214             : }
     215             : 
     216             : 
     217             : int
     218           0 : CharacteristicMap::getDomainDim() const {
     219           0 :     return domainDim;
     220             : }
     221             : 
     222             : 
     223             : int
     224           0 : CharacteristicMap::getImageDim() const {
     225           0 :     return imageDim;
     226             : }
     227             : 
     228             : 
     229             : std::vector<double>
     230        5280 : CharacteristicMap::eval(const std::vector<double>& ref_p, double eps) const {
     231        5280 :     if (static_cast<int>(ref_p.size()) != domainDim) {
     232           0 :         throw std::runtime_error("The argument's size doesn't match the domain"
     233           0 :                                  " dimension.");
     234             :     }
     235             : 
     236             :     // Find the nearest neighbor and its image values
     237             :     std::vector<int> nnIdxs;
     238        5280 :     if (findNearestNeighborIdxs(ref_p, nnIdxs, eps)) {
     239           0 :         return std::vector<double>(imageDim, std::stod("nan"));
     240             :     }
     241             :     // Image values of the nearest neighbor
     242        5280 :     const std::vector<double> y_nn = at(nnIdxs);
     243             :     // The result is based on the image values of the nearest neighbor
     244        5280 :     std::vector<double> y = y_nn;
     245             : 
     246             :     // Interpolate
     247       15840 :     for (int i = 0; i < domainDim; i++) {
     248             :         // Depending on the configuration of the points, different neighbors will be
     249             :         // used for interpolation
     250       10560 :         const double s = ref_p[i] - axes[i][nnIdxs[i]];
     251       10560 :         if (std::abs(s) <= eps) {
     252        6100 :             continue;
     253             :         }
     254        4460 :         bool b_constellation1 = s < 0 && nnIdxs[i] > 0;
     255             :         bool b_constellation2 = s >= 0
     256        2152 :                                 && nnIdxs[i] == static_cast<int>(axes[i].size()) - 1
     257        4460 :                                 && nnIdxs[i] > 0;
     258        2308 :         bool b_constellation3 = s < 0 && nnIdxs[i] == 0
     259        4460 :                                 && nnIdxs[i] < static_cast<int>(axes[i].size()) - 1;
     260             :         bool b_constellation4 = s >= 0
     261        4460 :                                 && nnIdxs[i] < static_cast<int>(axes[i].size()) - 1;
     262             : 
     263             :         double dx = 1;
     264             :         // Axis neighbor indices (i.e. the indices of the second support point)
     265        4460 :         std::vector<int> anIdxs = nnIdxs;
     266        4460 :         if (b_constellation1 || b_constellation2) {
     267        2308 :             anIdxs[i] -= 1;
     268        2308 :             dx = axes[i][nnIdxs[i]] - axes[i][anIdxs[i]];
     269        2152 :         } else if (b_constellation3 || b_constellation4) {
     270        2152 :             anIdxs[i] += 1;
     271        2152 :             dx = axes[i][anIdxs[i]] - axes[i][nnIdxs[i]];
     272             :         } else {
     273             :             continue;
     274             :         }
     275             :         // Image values of the axis neighbor
     276        4460 :         const std::vector<double> y_an = at(anIdxs);
     277             : 
     278        8920 :         for (int j = 0; j < imageDim; j++) {
     279             :             double dy = 0;
     280        4460 :             if (b_constellation1 || b_constellation2) {
     281        2308 :                 dy = y_nn[j] - y_an[j];
     282             :             } else {
     283        2152 :                 dy = y_an[j] - y_nn[j];
     284             :             }
     285             : 
     286             :             // Update
     287        4460 :             y[j] += s * dy / dx;
     288             :         }
     289             :     }
     290             : 
     291             :     return y;
     292             : }