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 319539 : CharacteristicMap::determineStrides() {
36 : strides.clear();
37 319539 : strides.reserve(domainDim);
38 319539 : strides.push_back(imageDim);
39 639078 : for (int i = 1; i < domainDim; i++) {
40 319539 : strides.push_back((int)axes[i - 1].size() * strides[i - 1]);
41 : }
42 319539 : }
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 : 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 : 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 : 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 : 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 : 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 319539 : CharacteristicMap::CharacteristicMap(const std::string& ref_mapString) {
142 : // Split the map string into its three main parts
143 958617 : const std::vector<std::string> tokens = StringTokenizer(ref_mapString, "|").getVector();
144 319539 : if (tokens.size() != 3) {
145 : 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 958617 : const std::vector<std::string> dimensionTokens = StringTokenizer(tokens[0], ",").getVector();
151 319539 : if (dimensionTokens.size() != 2) {
152 : throw std::runtime_error("The domain and image dimensions aren't specified"
153 0 : " correctly.");
154 : }
155 319539 : domainDim = std::stoi(dimensionTokens[0]);
156 319539 : imageDim = std::stoi(dimensionTokens[1]);
157 :
158 : // Create the map axes
159 958617 : const std::vector<std::string> axisTokens = StringTokenizer(tokens[1], ";").getVector();
160 319539 : if (static_cast<int>(axisTokens.size()) != domainDim) {
161 : throw std::runtime_error("The number of axes doesn't match the specified"
162 0 : " domain dimension.");
163 : }
164 958617 : for (auto& ref_axisToken : axisTokens) {
165 1917234 : std::vector<std::string> axisEntryTokens = StringTokenizer(ref_axisToken, ",").getVector();
166 : std::vector<double> axisEntries;
167 1919010 : for (auto& ref_axisEntryToken : axisEntryTokens) {
168 1279932 : axisEntries.push_back(std::stod(ref_axisEntryToken));
169 : }
170 639078 : axes.push_back(axisEntries);
171 639078 : }
172 :
173 : // Create the flattened map
174 958617 : const std::vector<std::string> flattenedMapTokens = StringTokenizer(tokens[2], ",").getVector();
175 319539 : int expectedEntryCnt = imageDim;
176 958617 : for (auto& ref_axis : axes) {
177 639078 : expectedEntryCnt *= (int)ref_axis.size();
178 : }
179 319539 : if (static_cast<int>(flattenedMapTokens.size()) != expectedEntryCnt) {
180 : 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 319539 : flattenedMap.reserve(expectedEntryCnt);
184 1627455 : for (auto& ref_flattenedMapToken : flattenedMapTokens) {
185 1307916 : flattenedMap.push_back(std::stod(ref_flattenedMapToken));
186 : }
187 :
188 319539 : determineStrides();
189 319539 : }
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 : 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 4460 : }
290 :
291 : return y;
292 5280 : }
|
