| 42 | } |
| 43 | |
| 44 | StoredFunction::SearchResult StoredFunction::search(double targetValue, double valueTolerance) const { |
| 45 | double minIndex = m_function.index(0); |
| 46 | double minValue = m_function.value(0); |
| 47 | |
| 48 | double maxIndex = m_function.index(m_function.size() - 1); |
| 49 | double maxValue = m_function.value(m_function.size() - 1); |
| 50 | |
| 51 | if (maxValue < minValue) { |
| 52 | std::swap(minIndex, maxIndex); |
| 53 | std::swap(minValue, maxValue); |
| 54 | } |
| 55 | |
| 56 | double index; |
| 57 | double value; |
| 58 | |
| 59 | if (targetValue < minValue) { |
| 60 | index = minIndex; |
| 61 | value = minValue; |
| 62 | } else if (targetValue > maxValue) { |
| 63 | index = maxIndex; |
| 64 | value = maxValue; |
| 65 | } else { |
| 66 | index = (minIndex + maxIndex) / 2; |
| 67 | value = m_function.interpolate(index); |
| 68 | |
| 69 | int searchDepth = 0; |
| 70 | |
| 71 | while ((std::fabs(targetValue - value) > valueTolerance) && (searchDepth < 64)) { |
| 72 | searchDepth++; |
| 73 | if (value < targetValue) { |
| 74 | minIndex = index; |
| 75 | minValue = value; |
| 76 | } else if (value > targetValue) { |
| 77 | maxIndex = index; |
| 78 | maxValue = value; |
| 79 | } |
| 80 | |
| 81 | double newIndex = (minIndex + maxIndex) / 2; |
| 82 | double newValue = m_function.interpolate(newIndex); |
| 83 | |
| 84 | // If at any point we move outside of the established upper and lower |
| 85 | // bound |
| 86 | // the function is not monotonic increasing or decreasing, and binary |
| 87 | // search |
| 88 | // can not be used so we have to bail out. |
| 89 | if (newValue > maxValue || newValue < minValue) |
| 90 | throw StarException("StoredFunction is not monotonic."); |
| 91 | |
| 92 | index = newIndex; |
| 93 | value = newValue; |
| 94 | } |
| 95 | } |
| 96 | |
| 97 | SearchResult result; |
| 98 | result.targetValue = targetValue; |
| 99 | result.searchTolerance = valueTolerance; |
| 100 | result.found = std::fabs(targetValue - value) <= valueTolerance; |
| 101 | result.solution = index; |