2
\$\begingroup\$

I'm trying to implement what I think griddedInterpolant(X,V,method) is doing inside Interp1.

I get the same results. Please review the extrapolation part as I'm not quite sure about it and comment about the algorithm correctness.

double Model::Interp1(const std::vector<double>& alpha, const std::vector<double>& varSurf, double detail_sensativity)
{
    //in case the interoplation point alraedy exists in the array
    //just return it
    std::vector<double>::const_iterator it = std::find(alpha.begin(), alpha.end(), detail_sensativity);
    if (it != alpha.end())
    {
        int index = it - alpha.begin();
        return varSurf[index];
    }
    else // we need linear interpolation/extrapolation.
    {
        double max = *std::max_element(alpha.begin(), alpha.end());
        double min = *std::min_element(alpha.begin(), alpha.end());
        if (detail_sensativity >= min && detail_sensativity <= max) //interpolate
        {
            for (size_t i = 0; i < alpha.size(); ++i)
            {
                if (detail_sensativity >= alpha[i] && detail_sensativity <= alpha[i + 1])
                {
                    //y = y0 + (y1-y0)*(x-x0)/(x1-x0);
                    double y0 = varSurf[i];
                    double y1 = varSurf[i + 1];
                    double x = detail_sensativity;
                    double x0 = alpha[i];
                    double x1 = alpha[i + 1];
                    double res = y0 + (y1 - y0)*(x - x0) / (x1 - x0);
                    return res;
                }
            }
        }
        else //extrapolate
        {
            //y = y0 + ((x - x0) / (x1 - x0)) * (y1 - y0)
            int i;
            if (detail_sensativity >= alpha.end)
            {
                i = varSurf.end() - varSurf.begin() - 1;                
            } 
            else
            {
                i = 0;
            }
            double y0 = varSurf[i];
            double y1 = varSurf[i + 1];
            double x = detail_sensativity;
            double x0 = alpha[i];
            double x1 = alpha[i + 1];
            double res = y0 + ((x - x0) / (x1 - x0)) * (y1 - y0);
            return res;
        }
    }
}
\$\endgroup\$
3
\$\begingroup\$

duplicated math

The two statements

double res = y0 + (y1 - y0)*(x - x0) / (x1 - x0);

and

double res = y0 + ((x - x0) / (x1 - x0)) * (y1 - y0);

are equivalent.

They both expand to

$$\frac{x y_0-x y_1+x_0 y_1-x_1 y_0}{x_0-x_1}$$

as can be seen in W\$\alpha \$ of the first and W\$\alpha \$ of the second (as the alternate forms are neatly in the same location, you can toggle between both tabs)

Even without the additional brackets around ((x - x0) / (x1 - x0)), it would still be the same thing.

That's somewhat plausible I'd say, because no matter if you are interpolating or extrapolating, you are still applying the formula for a line and that's the same whether you are between two points that define the line or outside of them.

duplicated code

For the code, that means the following block is a duplicate:

//y = y0 + (y1-y0)*(x-x0)/(x1-x0);
double y0 = varSurf[i];
double y1 = varSurf[i + 1];
double x = detail_sensativity;
double x0 = alpha[i];
double x1 = alpha[i + 1];
double res = y0 + (y1 - y0)*(x - x0) / (x1 - x0);
return res;

To simplify your code, move the block to the end. All code beforehand basically just figures out the right index i to use, then perform the code block. Here's an untested implementation of that. (only last else part of your code shown) // !!!!!!!!!! indicates lines of significant change.

else // we need linear interpolation/extrapolation.
{
    double max = *std::max_element(alpha.begin(), alpha.end());
    double min = *std::min_element(alpha.begin(), alpha.end());

    size_t i = 0;                                           // !!!!!!!!!!

    if (detail_sensativity >= min && detail_sensativity <= max) //interpolate
    {
        for (; i < alpha.size(); ++i)                       // !!!!!!!!!!
        {
            if (detail_sensativity >= alpha[i] && detail_sensativity <= alpha[i + 1])
            {
                break;                                      // !!!!!!!!!!
            }
        }
    }
    else //extrapolate
    {
        if (detail_sensativity >= alpha.end)
        {
            i = varSurf.end() - varSurf.begin() - 1;                
        }
    }

    double y0 = varSurf[i];
    double y1 = varSurf[i + 1];
    double x = detail_sensativity;
    double x0 = alpha[i];
    double x1 = alpha[i + 1];
    double res = y0 + ((x - x0) / (x1 - x0)) * (y1 - y0);
    return res;
}

Changes I made:

  1. moved intreaxtrapolation block to the end
  2. used size_t i to top and used it for both cases
  3. used break to end for loop. If you are a real C programmer, you'd shove all the condition of if (detail_sensativity >= alpha[i] && detail_sensativity <= alpha[i + 1]) into the condition of the for loop, just so you could have a loop without a body … just kidding, that wouldn't be very readable.

possible out of bounds error

Not sure if this is correct, just something looking suspicious:

Your loop ends at i < alpha.size();, thus i might end up being as big as alpha.size() - 1, which is the last valid index of the vector. However, you then do double x1 = alpha[i + 1];, which appears to be out of bounds in this case.

code organisation

You are writing a math utility function here. I would not put it into a Model namespace

         v--------------this
double Model::Inte....

I'm not a math expert, but I know that there are other ways to interpolate (cubic, bi-linear, Q-spline, slitherin, yoda, you name it). Maybe your Model asks for linear interpolation, but the function you are posting here is not a member function (assuming Model is a class). It doesn't depend on the state of an object and only uses the parameters passed to it (as far as I can tell). This looks like a free function and should possibly live in some Math namespace (or similar), independent from any Model that's using it.

Math is abstract and the same math is applied in many different fields in many different ways. Keep the abstraction by having the pure math functions separate and not tied to some Model that's interesting to your domain.

naming

It took me quite a while to figure out what the function is doing, because I was wondering what the heck detail_sensativity is supposed to be. Sure there's a typo in it, but still: what does it mean?

double Model::Interp1
(
    const std::vector<double>& alpha, 
    const std::vector<double>& varSurf, 
    double detail_sensativity
)

It's part of a bigger problem that's related to the former paragraph: you use domain specific names which are very cryptic to anybody outside your field.

I know what an interpolation function should be doing and with the clue that you are trying to duplicate griddedInterpolant(X,V,method) I knew what's going on. With your code alone, that would have been much harder. What is alpha? varSurf? detail_sensativity? And again, given there are many interpolation methods, what is Interp1?

Here's how this could look like:

double Interpolation::linear
(
    const std::vector<double>& x, 
    const std::vector<double>& y, 
    double x_value
)

simplification

I feel like this code could be simplified further. The case for extrapolation could be included in the for loop without being explicitly coded. After all, extrapolation just means using the same slope as the first and last interpolation steps. Checking max and min values feels like an overcomplication.

handling problematic input

At the moment, there are no error checks. One can happily pass in two std::vectors of different lengths and the algorithm will go bananas. There should be some error checks.

Also, it could make a lot more sense to use different input parameter types. There's a pair class for example. The x and y values should always come in pairs, so why not do that? (code untested, just to illustrate)

double Interpolation::linear
(
    const std::vector<std::pair<double, double>>& x_y_points,
    double x_value
)

Now you only have to check if the length is > 1. I mean you're having a Model there, so you probably model something. I think std::pair<double, double> models the fact that the values come in pairs pretty well. Definitely more so than having two vectors. You can always provide an overload for the naysayers.

\$\endgroup\$
2
\$\begingroup\$

Just a couple of minor additions to @I'll add comments tomorrow's excellent answer.

One simplification is to use std::minmax_element, to replace this:

double max = *std::max_element(alpha.begin(), alpha.end());
double min = *std::min_element(alpha.begin(), alpha.end());

with something like this:

auto mm = std::minmax_element(alpha.begin(), alpha.end());
min = *mm.first;
max = *mm.second;

This gives us a little extra efficiency by doing the search for both the minimum and maximum in a single sweep of the inputs. There might be some room for argument that the code is harder to read, but at least in my opinion not much--yes, it's one extra line of code, but at least if your input is large that's a pretty trivial price to pay for doubling the speed of this code.

Another simplification is in the code using std::find:

std::vector<double>::const_iterator it = std::find(alpha.begin(), alpha.end(), detail_sensativity);

This can be replaced with code like:

auto it = std::find(alpha.begin(), alpha.end(), detail_sensitivity);

This is one of the motivating examples for adding auto type deduction to C++ in the first place. Unless you're using a really old compiler, you might as well take advantage of it.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.