Fill dynamic size vector recursively

Question

Maybe let me state my situation in pseudo-C++-code first:

std:vector<double> sample(someFunctor f, double lower, double upper) {
    double t = (lower + upper)/2;
    double newval = f(t);

    if (f(upper) - newval > epsilon)
        subsample1 = sample(f, t, upper);
    if (newval - f(lower) > epsilon)
        subsample2 = sample(f, lower, t);

    return concat(subsample2, newval, subsample1);
}

where concat just, well, concats the returned vectors. Basically, I am sampling a function in a way such that between two saved function values there is only a small difference.

I am not happy with the way stated above as there seems to be quite a bit of memory allocation (allocate two subvectors, then concat those and another element) in every recursive step. That piece of code will have to run in a part of my algorithm that is critical with respect to performance. Once upper - lower is fairly small, evaluating f will not take a great amount of time.

So my questions:

• Do you see a clever way to use the same data structure in all recursive calls and just fill the current parts of that vector? (Keeping in mind that the number of function evaluations needed is not known upfront). Thoughts on this:

  • Using a list instead of a vector. But I feel the memory overhaul is not adequate for just storing doubles.
  • Keep holes in the vector and maintain another vector saying which entries are filled. A the end of a recursive call shift the entries so that there are no holes between the subsamples and newval. But now I switch copying by shifting with additional work for the second vector - probably bad idea.

• Do you see a way to get rid of the recursion totally? However, for correctness it is important that I use the divide-and-conquer pattern stated above. The function f makes heavy use of the upper and lower bounds and gains quite a big deal of performance by this.

Thank's for your thoughts.

---

As per Space_C0wb0y's request let me try to rephrase my problem. Maybe the first explanation was not very clear.

I have some function (in a mathematical sense) that I want to sample (e.g. to evaluate at certain points) in a given interval.

Suppose that interval is [0,100]. I know the functions values at 0 and at 100. Maybe that is f(0)=0 and f(100) = 40. Now I evaluate the function at the intervals midpoint, which is 50. Say, my function returns f(50)=10. As f(0)-f(50) <= 10, I do not need further samples in the interval [0,50]. However, I need further computations for the interval [50,100]. Thus, in the next (recursive) step I evaluate f(75). Now repeat the above logic recursively.

At the end I would like to (two) vectors giving me the function values with the corresponding parameters like this:

parameter  = vector(0, 50, 56.25, 62.5, 75, 100)
value      = vector(0, 10, 17.21, 25    34,  40)

I am looking for the best (as is most performant) approach to build these vectors recursively.

Hope this clarifies things.

Answer 1

Since space is not your major concern so I'll go on using the recursion.

1. Use copy by reference instead of copy by (return) value.

2. No need to pass in functor as it's constant.

3. It could have been faster if low and high are integers. It's depends on requirements though.

    // Thanks to Space_C0wb0y, here we avoid using a global vector
    // by passing the vector as reference. It's efficient as there
    // is no copy overhead as well.        
    void sample(vector<double>& samples, double low, double high)
    {
       // You can use shift operator if they're integers.
       double mid = (low + high)/2;

       // Since they're double, you need prevent them from being too close.
       // Otherwise, you'll probably see stack overflow.
       // Consider this case:
       // f(x): x=1, 0<x<8;  x*x, x<=0 or x>=8
       // low = 1, high = 10, epsilon = 10
       if (high - low < 0.5)
       {
          samples.push_back(f(mid));
          return;
       }   

       // The order you write the recursive calls guarantees you
       // the sampling order is from left to right.
       if (f(mid) - f(low) > epsilon)
       {
          sample(samples, low, mid);
       }

       samples.push_back(f(mid));

       if (f(high) - f(mid) > epsilon)
       {
          sample(samples, mid, high);
       }   
    }

Answer 2

I would recommend the following approach:

1. Don't use two vectors, but rather one vector with pairs or a custom struct to represent parameters and values:

   struct eval_point {
       double parameter;
       double value;
   };
   
   std::vector<eval_point> evaluated_points;
   

2. Change your algorithm to write the results of the evaluations to an output iterator:

   template<class F, class output_iterator_type>
   void sample(F someFunctor, double lower, double upper,
               output_iterator_type out) {
       double t = (lower + upper)/2;
       eval_point point = { t, f(t) };
   
       if (f(upper) - point.value > epsilon) {
           *out = point;
           ++out;
           sample(f, t, upper, out);
       }
       if (point.value - f(lower) > epsilon) {
           *out = point;
           ++out;
           subsample2 = sample(f, lower, t, out);
       }
   }
   

   The above is a modification of your pseudo-code showing what it might look like when using an output iterator. It is not tested, so I am not sure if it is correct. In principle, you would call it like this:

   std::vector<eval_point> results;
   sample(someFunction, 0, 100, std::back_inserter<eval_point>(results));
   

   This way you will not have to create new vectors for each recursive call. If you can guess a reasonable lower bound for the number of samples, you might be able to preallocate such that no reallocations are required. In that case you would call it like this:

   std::vector<eval_point> results(lower_bound_for_samples);
   sample(someFunction, 0, 100, results.begin());
   

You would then have to add an extra counter to keep track of how many samples were generated.

Answer 3

I don't see why you decline a list solution. The worst case would be that your list has 3 times the size than your raw data. I think that is by far less than what you have when you create a new vector on each function-call. You should try it out as it does not require that much change, because the interface of both is nearly the same.