Skip to content

Adding missing examples for Tree Traversal and FFT in C++. #114

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed

Conversation

Gathros
Copy link
Contributor

@Gathros Gathros commented May 16, 2018

I added DFT to FFT, DFS_recursive_postorder and DFS_recursive_inorder_btree to tree traversal all in C++.

@strega-nil
Copy link

strega-nil commented May 16, 2018

I don't like this change. It looks very C-y, and there is bit operations on int, which doesn't seem right. What do we actually have to change?

@Gathros
Copy link
Contributor Author

Gathros commented May 16, 2018

We need to add dft in fft.cpp, dfs_recursive_postorder and dfs_recursive_inorder_btree in tree_example.cpp.

@strega-nil strega-nil self-requested a review May 16, 2018 15:48
@strega-nil
Copy link

strega-nil commented May 16, 2018

I'm going to add a PR to your PR. Thanks for building the base!

@Gathros Gathros added the Implementation Edit This provides an edit to an algorithm implementation. (Code and maybe md files are edited.) label May 16, 2018
@strega-nil
Copy link

@Gathros poke?

@zsparal
Copy link
Contributor

zsparal commented May 20, 2018

I'd argue that the use of iterators and std::accumulate actually reduces the readability of the code. It would be fine if we wanted to make a standard library implementation but it's not very suitable for an educational resource. Also, having to capture mutable variables in an accumulate is a pretty good indication you shouldn't be using it in my opinion. I think this version is a lot more readable:

std::vector<c64> discrete(const std::vector<c64>& vec) {
    std::vector<c64> result(vec.size());

    auto size = vec.size();
    for (size_t i = 0; i < size; ++i) {
        for (size_t j =0 ; j < size; ++j) {
            result[i] += vec[j] * std::exp(c64(0, -2.0 * pi<double>() * j * i / size));
        }
    }

    return result;
}

@zsparal
Copy link
Contributor

zsparal commented May 20, 2018

Also, I'll say it again: std::valarray is literally made for vector math like this, so if we want to be really gung-ho about using standard library features everywhere then using it would make the code a lot shorter and more readable (than reimplementing all their features with iterators/std::vectors)

auto v = c64(1.0);
for (size_t k = 0; k < stride / 2; k++) {
for (std::size_t k = 0; k < stride / 2; k++) {
first[k + j + stride / 2] =
first[k + j] - v * first[k + j + stride / 2];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This indentation looks incorrect

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what is incorrect about it?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh sorry, I actually misread that line

@strega-nil
Copy link

strega-nil commented May 20, 2018

@Gathros I'm not sure I agree. I prefer the one with the standard algorithms, although I think some of it could be worked around better (i.e., give names)

@strega-nil
Copy link

Also, I just realized what that call to exp is doing, and that seems like a really badly worded way of doing a sine or cosine

@Gathros
Copy link
Contributor Author

Gathros commented May 20, 2018

It's @Gustorn not me @ubsan.

@strega-nil
Copy link

Sorry

void cooley_tukey(Iter first, Iter last) {
auto size = last - first;
if (size >= 2) {
void recursive(Iter first, Iter last) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As I said in another comment, I'd still argue for valarray here, it makes the implementation a lot more readable:

void fft(std::valarray<c64>& v) {
    auto size = v.size();
    if (size < 2) {
        return;
    }

    auto even = std::valarray<c64>(v[std::slice(0, size / 2, 2)]);
    auto odd  = std::valarray<c64>(v[std::slice(1, size / 2, 2)]);

    fft(even);
    fft(odd);

    for (size_t k = 0; k < size / 2; k++) {
        auto w = std::polar(1.0, -2.0 * pi<double> * k / size) * odd[k];
        v[k] = even[k] + w;
        v[k + size / 2] = even[k] - w;
    }
}

Copy link

@strega-nil strega-nil May 20, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems like that's more of an issue with the current implementation; we could do

template <typename Iter>
void recursive(Iter first, Iter last) {
  std::size_t size = last - first;
  if (size < 2) {
    return;
  } else if (size % 2 != 0) {
    std::cerr << "The range must be a power of two size\n";
    std::abort();
  }

  auto odds = std::vector<c64>(size / 2);
  auto evens = std::vector<c64>(size / 2);

  for (std::size_t i = 0; i < size / 2; ++i) {
    odds[i] = first[i * 2 + 1];
    evens[i] = first[i * 2];
  }

  // recurse the splits and butterflies in each half of the range
  recursive(begin(odds), end(odds));
  recursive(begin(evens), end(evens));

  // now combine each of those halves with the butterflies
  for (std::size_t k = 0; k < size / 2; ++k) {
    auto w = std::polar(1.0, -2.0 * pi<double>() * k / size) * odds[k];
    first[k] = evens[k] + w;
    first[k + size / 2] = evens[k] - w;
  }
}

// versus

void recursive(std::valarray<c64>& v) {
  auto size = v.size();
  if (size < 2) {
    return;
  } else if (size % 2 != 0) {
    std::cerr << "The range must be a power of two size\n";
    std::abort();
  }

  auto even = std::valarray<c64>(v[std::slice(0, size / 2, 2)]);
  auto odd = std::valarray<c64>(v[std::slice(1, size / 2, 2)]);

  recursive(even);
  recursive(odd);

  for (size_t k = 0; k < size / 2; k++) {
    auto w = std::polar(1.0, -2.0 * pi<double>() * k / size) * odd[k];
    v[k] = even[k] + w;
    v[k + size / 2] = even[k] - w;
  }
}

which is only 4 lines of useful extra code; I don't think the latter is significantly more readable, and I actually prefer oddsₖ = first₂ₖ₊₁ , evensₖ = first₂ₖ

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is even nicer; moving out the split into it's own function, to make it really clear what's happening:

template <typename In, typename Out0, typename Out1>
void split(In initial, In const last, Out0 out0, Out1 out1) {
  if ((last - initial) % 2 != 0) {
    std::cerr << "The range must be a power of two size\n";
    std::abort();
  }

  while (initial != last) {
    *out0 = *initial;
    ++out0;
    ++initial;
    *out1 = *initial;
    ++out1;
    ++initial;
  }
}

// `recursive` does the cooley-tukey algorithm, recursively
template <typename Iter>
void recursive(Iter const first, Iter const last) {
  std::size_t size = last - first;
  if (size < 2) {
    return;
  }

  auto evens = std::vector<c64>(size / 2);
  auto odds = std::vector<c64>(size / 2);
  split(first, last, begin(evens), begin(odds));

  // recurse the splits and butterflies in each half of the range
  recursive(begin(odds), end(odds));
  recursive(begin(evens), end(evens));

  // now combine each of those halves with the butterflies
  for (std::size_t k = 0; k < size / 2; ++k) {
    auto w = std::polar(1.0, -2.0 * pi<double>() * k / size) * odds[k];
    first[k] = evens[k] + w;
    first[k + size / 2] = evens[k] - w;
  }
}

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I disagree with the separate split version being nicer, see my comment in the main thread. We're not really writing a standard library grade implementation and split just obfuscates a very significant step in the algorithm. I do like the proposed odds/evens version however, although I would still go with the iterator-free version with a proper return value:

void fft_calculate_recursive(std::vector<c64>& values) {
    std::size_t size = values.size();
    if (size < 2) {
        return;
    } else if (size % 2 != 0) {
        std::cerr << "The range must be a power of two size\n";
        std::abort();
    }

    auto odds = std::vector<c64>(size / 2);
    auto evens = std::vector<c64>(size / 2);

    for (std::size_t i = 0; i < size / 2; ++i) {
        odds[i] = values[i * 2 + 1];
        evens[i] = values[i * 2];
    }

    // recurse the splits and butterflies in each half of the range
    fft_calculate_recursive(odds);
    fft_calculate_recursive(evens);

    // now combine each of those halves with the butterflies
    for (std::size_t k = 0; k < size / 2; ++k) {
        auto w = std::polar(1.0, -2.0 * pi<double>() * k / size) * odds[k];
        values[k] = evens[k] + w;
        values[k + size / 2] = evens[k] - w;
    }
}

std::vector<c64> recursive(const std::vector<c64>& values) {
  std::vector<c64> result(values.begin(), values.end());
  fft_calculate_recursive(result);
  return result;
}

I think that teaching people to return actual results rather than using out parameters is way more important than any perceived genericity one gets from iterators (I'm not arguing that iterators aren't more general - I'm arguing that it isn't needed here).

@zsparal
Copy link
Contributor

zsparal commented May 20, 2018

I read modern C++ on a regular basis and I needed to take a second look at that generate/aggregate combo to realize what it's doing. aggregate is not the algorithm to use when you need external state. Not only does it make it much harder to grasp what the algorithm is actually doing, it also ruins thread safety (if you ever wanted to switch it to a reduce for example). With range-v3 we could write this nicely but with the existing STL algorithms? Not so much.

@strega-nil
Copy link

@Gustorn I have figured out how to do it nicely, one sec

@strega-nil
Copy link

strega-nil commented May 20, 2018

@Gustorn How's this?

template <typename F>
auto sum(std::size_t const start, std::size_t const end, F f) {
  using return_type = decltype(f(start));
  auto ret = return_type(0);

  for (auto i = start; i < end; ++i) {
    ret += f(i);
  }

  return ret;
}

template <typename It, typename F>
It indexed_generate(
    It first, std::size_t const start, std::size_t const end, F f) {
  auto it = first;
  for (auto i = start; i < end; ++i, ++it) {
    *it = f(i);
  }
  return it;
}

template <typename Iter>
void discrete(Iter first, Iter last) {
  auto size = last - first;
  auto temp = std::vector<c64>();

  indexed_generate(std::back_inserter(temp), 0, size, [&](auto i) {
    return sum(0, size, [&](auto k) {
      // NOTE:
      // this is equivalent to cos(-2 pi i k / size) + i sin(-2 pi i k / size)
      return first[k] * std::exp(c64(0, -2.0 * pi<double>() * i * k / size));
    });
  });

  std::copy(begin(temp), end(temp), first);
}

(I'm not sure what to call indexed_generate, but you get the idea)

We want to represent the idea that temp is defined by the function:

aᵢ = Σ(k = 0 to N) bₖ (cos(-2πkj / N) + i sin(-2πkj / N))

where b is the input function.

@strega-nil
Copy link

And a second attempt:

template <typename Iter>
void discrete(Iter first, Iter last) {
  using namespace std::literals;

  auto const original = std::vector<c64>(first, last);
  auto const size = original.size();

  auto const i2pi = -2.0i * pi<double>();

  for (std::size_t i = 0; i < size; ++i) {
    first[i] = sum(0, size, [&](auto k) {
      return original[k] * std::exp(i2pi * double(i * k) / double(size));
    });
  }
}

@zsparal
Copy link
Contributor

zsparal commented May 20, 2018

Yeah, I like the last version a lot more. You could also do something like:

template <typename Iter, typename F>
auto sum(Iter first, Iter last, const F& fn) {
    using namespace std::literals;
    return std::accumulate(first, last, std::make_pair(0, 0i), [&fn](auto acc, const auto& x) {
        return std::make_pair(acc.first + 1, acc.second + fn(acc.first, x));
    }).second;
}

template <typename Iter>
void discrete(Iter first, Iter last) {
  using namespace std::literals;

  auto const original = std::vector<c64>(first, last);
  auto const size = original.size();

  for (std::size_t i = 0; i < size; ++i) {
    first[i] = sum(original.begin(), original.end(), [i, size](auto k, auto value) {
      return value * std::exp(-2.0i * pi<double>() * double(i * k) / double(size));
    });
  }
}

Alternatively for sum (I actually prefer this one):

template <typename Iter, typename F>
auto sum(Iter first, Iter last, const F& fn) {
    auto result = c64(0, 0);
    size_t i = 0;
    for (; first != last; ++first) {
        result += fn(i++, *first);
    }
    return result;
}

@strega-nil
Copy link

strega-nil commented May 20, 2018

Alright, then we can change over to

template <typename It, typename F>
auto sum(It const first, It const last, F f) {
  using return_type = decltype(f(0, *first));
  auto ret = return_type(0);

  std::size_t i = 0;
  for (auto it = first; it < last; ++it, ++i) {
    ret += f(i, *it);
  }

  return ret;
}

template <typename Iter>
void discrete(Iter const first, Iter const last) {
  using namespace std::literals;

  auto const original = std::vector<c64>(first, last);
  auto const size = original.size();

  auto const i2pi = -2.0i * pi<double>();

  for (std::size_t i = 0; i < size; ++i) {
  for (std::size_t i = 0; i < size; ++i) {
    first[i] =
        sum(begin(original), end(original), [&](auto const k, auto const el) {
          return el * std::exp(i2pi * double(i * k) / double(size));
        });
  }

  }
}

@Gathros

@Gustorn I still believe that we should be teaching iterators - the great thing about them is that we can switch out std::array<c64, 64> for std::valarray<c64> without changing any other code shrug (also, your implementation doesn't work)

@Gathros
Copy link
Contributor Author

Gathros commented May 20, 2018

@ubsan Changed it.

@zsparal
Copy link
Contributor

zsparal commented May 20, 2018

Sorry, I only checked if it compiles, not if it works (on godbolt, didn't have a C++ compiler handy). The point I was making with valarray is that you can write the algorithm in an easier and more understandable manner. Yes, the generic iterator version will work but the valarray one is nicer.

As for iterators: I'm not the biggest fan of them in general. If we want to be generic, I'd rather be generic over the container and not the iterators. That would still allow you to use range-based for loops and standard algorithms (and that's kind of what the standard is moving towards with ranges anyway).

@zsparal
Copy link
Contributor

zsparal commented May 21, 2018

I'll try to make one last argument against the iterator/algorithm solution. First, let's compare the now final version with my proposed one:

template <typename It, typename F>
auto sum(It const first, It const last, F f) {
  using return_type = decltype(f(0, *first));
  auto ret = return_type(0);

  std::size_t i = 0;
  for (auto it = first; it < last; ++it, ++i) {
    ret += f(i, *it);
  }

  return ret;
}

template <typename Iter>
void discrete(Iter const first, Iter const last) {
  using namespace std::literals;

  auto const original = std::vector<c64>(first, last);
  auto const size = original.size();

  auto const i2pi = -2.0i * pi<double>();
  for (std::size_t i = 0; i < size; ++i) {
    first[i] =
        sum(begin(original), end(original), [&](auto const k, auto const el) {
          return el * std::exp(i2pi * double(i * k) / double(size));
        });
  }

  }
}

vs

std::vector<c64> discrete(const std::vector<c64>& vec) {
    std::vector<c64> result(vec.size());

    auto size = vec.size();
    for (size_t i = 0; i < size; ++i) {
        for (size_t j =0 ; j < size; ++j) {
            result[i] += vec[j] * std::exp(c64(0, -2.0 * pi<double>() * j * i / size));
        }
    }

    return result;
}

I don't know about you, but I find the latter version immensely more readable. I understand the notion of avoiding raw loops but in this case there's one line of self-contained math code. There are some problems that benefit from abstractions, this isn't one of them. The proposed implementation also doesn't really give that much genericity over my proposed solution (it requires a random access iterator). If we really want to support multiple containers we could just do:

template <typename T>
T discrete(const std::vector<c64>& vec) {
    T result(vec.size());

    auto size = vec.size();
    for (size_t i = 0; i < size; ++i) {
        for (size_t j =0 ; j < size; ++j) {
            result[i] += vec[j] * std::exp(c64(0, -2.0 * pi<double>() * j * i / size));
        }
    }

    return result;
}

The problem with always going with the "modern C++" solution without trying to evaluate to potential benefits is that it introduces a lot of noise, most of which is completely unnecessary for the algorithm at hand.

Another important thing is that this book is supposed to teach algorithms, not C++. That being said, I still think we should write clean, readable and performant C++ code but we shouldn't necessarily strive for standard library-like implementations. Most people will never need the genericity it provides and it obfuscates the important parts of the algorithms. Also, chances are that people who do need that level of genericity can easily adapt the container version to suit their needs. The same cannot be said for beginners who're reading about an algorithm for the first time.

@strega-nil
Copy link

strega-nil commented May 21, 2018

@Gustorn In fact, I disagree completely. The latter two are not more readable, and this book is attempting to teach good C++. If you don't want to read C++, you have the choice to use any of the other languages. We are attempting to represent the idea

capture

in code, and so we want to have the closest to that idea possible, preferably using algorithms.

I could see the argument for using a specific data type, instead of iterators, but stuff like sum is so much more readable and understandable than the for-loop with a plus-equals.

My opinion is that writing code like that gives people a surface level understanding on what's going on, without giving any deeper mathematical meaning to the code - it's clear, for example, what bit_reverse_sort is doing mechanically, yet I have zero idea what it's physically doing (which is why it needs comments).

@zsparal
Copy link
Contributor

zsparal commented May 21, 2018

Alright, I was trying to be a little more subtle about it but the new code isn't good C++. It's C++ that might look good in a CppCon talk that showcases language features.

Sure, you might've given a nice name to sum but it doesn't even just sums the elements.. It's basically a map and a sum and with the use of auto in the lambda parameters you don't even know which parameter is supposed to be the index.

On the other hand you have a for loop with a single state changing operation which is +=. If that doesn't scream sum to you I don't know what to tell you. The thing you seem to be forgetting is that sum is part of the algorithm so the reader will have to understand that as well. Let's take a look at the first few lines of the function:

  using return_type = decltype(f(0, *first));
  auto ret = return_type(0);

You very rarely have to use this pattern in everyday code. This is the kind of code that gives modern C++ advocates a bad name (which is a shame because I'm one of them when it makes sense): just add genericity and abstractions to a problem until it becomes completely unrecognizable. Do we want to emphasize that it's a sum? Sure, we can do that with some straightforward code:

c64 fft_sum(size_t i, const std::vector<c64>& values) {
    auto sum = c64(0);
    auto size = values.size();
    for (size_t k =0 ; k < size; ++k) {
        sum += values[k] * std::exp(c64(0, -2.0 * pi<double>() * k * i / size));
    }
    return sum;
}

std::vector<c64> discrete(const std::vector<c64>& values) {
    std::vector<c64> result(values.size());

    auto size = values.size();
    for (size_t i = 0; i < size; ++i) {
        result[i] = fft_sum(i, values);
    }

    return result;
}

I'm still against this version because it still obfuscates part of the algorithm but at least it doesn't add unneeded complexity.

@strega-nil
Copy link

strega-nil commented May 21, 2018

@Gustorn I disagree absolutely completely. We should write what we want to happen, in a declarative style, not in an imperative, tell the machine what should happen, sort of way. This is the basis of good code design.

This is why I prefer the original sum - I want to point out that this is a sigma operation. I also frankly would prefer it if we used c64(cos(-2πkj), sin(-2πkj)), because that actually shows what we mean.

In a perfect world, I would write:

template <typename It>
std::vector<c64> dft(It const original, It const last) {
  std::size_t const size = last - original;
  auto ret = std::vector<c64>(size);

  std::iota_transform(begin(ret), end(ret), [original] (std::size_t n) {
    return sum(0, size, [original] (std::size_t k) {
      return original[k] * c64(
          std::cos(-2.0 * pi * n * k / size),
          std::sin(-2.0 * pi * n * k / size));
    }
  });

  return ret;
}

Which is basically an exact translation of the mathematical operation.

The operation is not mutating the values of a vector, and returning it. The operation is a function from the naturals to the complex reals, where each value is a sum. You are obfuscating the true operation by introducing mutation where it should be declarative.

@zsparal
Copy link
Contributor

zsparal commented May 21, 2018

Yup, the pinnacle of declarative code right here:

template <typename It, typename F>
auto sum(It const first, It const last, F f) {
  using return_type = decltype(f(0, *first));
  auto ret = return_type(0);

  std::size_t i = 0;
  for (auto it = first; it < last; ++it, ++i) {
    ret += f(i, *it);
  }

  return ret;
}

As I mentioned, sum is part of the algorithm. This is why abstraction doesn't work well here: you have to understand the whole code to understand how it works, you're not just going to look at pieces of it in isolation. Also the whole notion that declarative style is always better than imperative style is kind of ridiculous. The key is to recognize which one to use and when. All that being said, I'm not very motivated to keep going with this argument since it feels like I'm arguing against C++ blog posts and talks and not against actual technical arguments.

@strega-nil
Copy link

strega-nil commented May 21, 2018

@Gustorn the point is that sum is an implementation function. That's why it has a name. You don't need to read it. The fact that it's not in the standard library, frankly, is kind of silly.

Also, for clearly mathematical operations like this, of course declarative is better than imperative. Imperative code just hides the underlying algorithm

@strega-nil
Copy link

strega-nil commented May 21, 2018

Even if we don't end up with iterator functions, we can write:

std::valarray<c64> recursive(std::valarray<c64> input) {
  auto const size = input.size();
  if (size < 2) {
    return input;
  }

  auto evens = std::valarray<c64>(input[std::slice(0, size / 2, 2)]);
  auto odds = std::valarray<c64>(input[std::slice(1, size / 2, 2)]);

  evens = recursive(evens);
  odds = recursive(odds);

  // now combine each of those halves with the butterflies
  auto ret = std::valarray<c64>(size);
  for (std::size_t k = 0; k < size / 2; ++k) {
    auto w = odds[k] * std::polar(1.0, -2.0 * pi<double>() * k / size);
    ret[k] = evens[k] + w;
    ret[k + size / 2] = evens[k] - w;
  }
  return ret;
}

// `iterative` does the cooley-tukey algorithm iteratively
std::valarray<c64> iterative(std::valarray<c64> input) {
  sort_by_bit_reverse(begin(input), end(input));

  // perform the butterfly on the range
  auto const size = input.size();
  for (std::size_t stride = 2; stride <= size; stride *= 2) {
    auto w = exp(c64(0, -2.0 * pi<double>() / stride));
    for (std::size_t j = 0; j < size; j += stride) {
      auto v = c64(1.0);
      for (std::size_t k = 0; k < stride / 2; k++) {
        input[k + j + stride / 2] =
            input[k + j] - v * input[k + j + stride / 2];
        input[k + j] -= (input[k + j + stride / 2] - input[k + j]);
        v *= w;
      }
    }
  }

  return input;
}

std::valarray<c64> discrete(std::valarray<c64> input) {
  auto const size = input.size();

  auto ret = std::valarray<c64>(size);

  iota_transform_n(begin(ret), 0, size, [&](std::size_t i) {
    return sum(0, size, [&](auto const k) {
      auto const coefficient = input[k];
      auto const expression =
          c64(std::cos(-2.0 * pi<double>() * i * k / size),
              std::sin(-2.0 * pi<double>() * i * k / size));
      return coefficient * expression;
    });
  });

  return ret;
}

which is much more understandable, because it actually declaratively says "this is the math that we are doing".

@leios
Copy link
Member

leios commented Jun 28, 2018

Hey @ubsan, I am sorry it has taken me so long to get to this PR. This is something we have discussed in detail internally on discord and through a few other platforms, but we feel that the code submitted to the algorithm archive in C++ is incredibly difficult to read and follow. This is mostly due to the syntactical complexity of the code, which does not seem standard, even for well-seasonsed C++ programmers.

As I've said before, you know your stuff, that is clear; however, if individuals who have been reading and writing C++ code for years have trouble understanding the code here, then there is a problem. At first, I thought it was just my own misunderstandings because my C++ code is basically C; however, it has come to my attention that other members of the community agree that the code is too difficult to read, and this is a problem.

I want to make clear that we definitely appreciate the help and support and would love to see you stay a member of the community and continue to review and write code for the AAA (we are not kicking you out); however, tone down the syntactical complexity to match what we have said before in the style guideline issue (#18) and How To Contribute section, which basically say that the code written should be "should be readable and understandable to anyone -- especially those who are new to the language." Your code definitely does not fit this criteria.

The code provided in #18 was definitely clearer than code provided here, and that should be roughly the syntactical complexity of future code in the archive.

@strega-nil
Copy link

@leios I disagree completely. If you want imperative code, you can write it yourself. The code as you've written it in Julia utterly obsfucates the underlying mathematical concepts; if I'm going to write code for something like this, I want to be able to see the math underlying it - I don't want to see the mutations. I'm unwilling to write the code in an imperative manner.

@leios
Copy link
Member

leios commented Jun 29, 2018

You are right, I have a bias towards imperative code; however, @ubsan 3 things:

  1. I should have been more clear: this isn't about the FFT, specifically. It was a general note that the code you have written for the AAA is incredibly over-engineered every time. There have been a few cases where I genuinely needed an implementation in C++ for certain algorithms and I could not use the implementation in the Algorithm Archive for C++ because it was unreadable and unusable in my codebase. I instead just used the C code with small modifications. That should not have happened.

I know what you are going to say, "Well, any good C++ programmer should be able to understand this code," and you are not necessarily wrong. The problem is that the people who need the Algorithm Archive are people like me -- those who use C++ frequently enough that they need an algorithm implemented in the language, but are not "good" enough to take the mathematical definition and immediately write it in code. They need code they can quickly look at and say, "Ah. That's how you do it. I bet I can get this working this afternoon!"

  1. I am not sure how the code I wrote "utterly obfuscates the underlying mathematical concepts," simply because it is written in an imperative style. As mentioned before (by Gustorn ~May 21st), there will always be slight variations of the mathematical definition and algorithmic implementation. If you want math, you read the chapter. If you want a language implementation, you read the code. Some people (like me) understand better when I see a few for loops, because it helps them grasp the dimensionality of the problem better than mathematical notation. If you want to write in a different coding style, pick a language that better supports that, or find ways to make it syntactically clearer what you are actually doing.

This is (decently) clear:

template <typename T>
T discrete(const std::vector<c64>& vec) {
    T result(vec.size());

    auto size = vec.size();
    for (size_t i = 0; i < size; ++i) {
        for (size_t j =0 ; j < size; ++j) {
            result[i] += vec[j] * std::exp(c64(0, -2.0 * pi<double>() * j * i / size));
        }
    }

    return result;
}

Sure, people who don't know C++ probably don't like the template, but that's fine. The std::exp(c64(0, -2.0 * pi<double>() * j * i / size)) is also a bit weird, but that's good ol' C++ for you.

This is a short sample of some of the code you wrote:

template <typename It>
std::vector<c64> dft(It const original, It const last) {
  std::size_t const size = last - original;
  auto ret = std::vector<c64>(size);

  std::iota_transform(begin(ret), end(ret), [original] (std::size_t n) {
    return sum(0, size, [original] (std::size_t k) {
      return original[k] * c64(
          std::cos(-2.0 * pi * n * k / size),
          std::sin(-2.0 * pi * n * k / size));
    }
  });

  return ret;
}

I can get a few lines in... Again, templates are weird, but it is what it is. Your constant use of auto makes auto ret = std::vector<c64>(size); a bit odd to read (why not just std::vector<c64> ret(size)?), but no big deal there. My main questions come after this. What is std::iota_transform(...)? Google doesn't know either. Is that a lambda somewhere hidden in the middle? Where is sum defined? Wait, is that a lambda in a lambda?

As far as using Euler's formula, which you seemed to dislike... Well, it's how the FT and FFT are defined. It's actually a bit misleading in my opinion to use sines and cosines unless your language demands it. Again, though, this is kinda personal preference. Everyone knows the formula, so either will do.

I don't mind you writing code in a different style. C++ supports a number of different styles. It's just that your code is so syntactically complex that style doesn't matter. It is simply hard to read and understand.

Disclaimer: I am aware that these are just snippets from the above conversation. They were not chosen delicately to prove my point or anything... they were just the smallest chunks of code I could grab that both seemed to be doing the same thing.

  1. You are right, I should re-write the julia code. It's a bit messy in certain places. I'll see about fixing that soon. Thanks for bringing that to my attention!

@strega-nil
Copy link

strega-nil commented Jun 29, 2018

@leios It's really not complex. It's not what you're used to, maybe, but it's fairly close to what I'd write in, say, OCaml.

val dft : Complex.t array -> Complex.t array

let dft original =
  let size = Array.length original in
  Array.init size (fun n ->
    sum 0 size (fun k ->
      let c = -2.0 *. pi *. (float_of_int n) *. (float_of_int k) /. (float_of_int size) in
      Complex.(original.(k) * { re= cos c; im= sin c }) ))

It's how one writes it in a non-imperative style, at least imo.

@jiegillet
Copy link
Member

If I may butt in, I do like that piece of OCaml, it's quite readable, but the C++ version is not readable at all. C++ is not very suited to declarative style, you have to jump through hoops to write the code you like. Maybe submit some OCaml code instead? ^^

@leios
Copy link
Member

leios commented Jun 29, 2018

Yeah, I was just about to write a similar comment. That code seems fine to me.

@strega-nil
Copy link

strega-nil commented Jun 29, 2018

They're... exactly the same code. What?

Literally the only difference is that ret gets declared first in C++, and then written into with iota_transform (a function I invented as an example), rather than created in the call to Array.init.

@jiegillet
Copy link
Member

They do exactly the same thing, but the OCaml is more readable, because the syntax is natural. For example, lambda functions. I don't know OCaml at all, and my C++ is worthless, but reading your OCaml, I see straight away that the syntax is fun x -> x. I've read your C++ 12 times and I can't figure out the syntax. It's just hard to read and not a natural feature of the language.

@strega-nil
Copy link

I don't know how to argue with that. I think I'm done with this thread. It seems to me we'll get yet another C-with-classes resource online. yay

@leios
Copy link
Member

leios commented Jun 29, 2018

@ubsan Thanks for answering. I think I have a better understanding of where you are coming from and the general differences in style between you and other coders here.

Again, to be clear: I don't think you are wrong. I think that C++ allows for the style you want, and that's fine. To be honest, I have genuinely learned a lot from reading through your code.

My understanding is that you would like to stay away from imperative code so you can be closer to the mathematical foundations of the algorithms you are using. That makes a lot of sense; however, we are arguing that in doing so, the C++ code becomes hard to read. This is probably because most people still use C++ imperatively, even if it is "multi-paradigmatic."

I feel there must be common ground. We are not doing completely OO code here, either. We use what makes sense, and the community dictates what makes sense in each case.

We have had a bit more discussion on this thread, so let's see if anyone else in the community chimes in. Maybe we should create a separate issue?

dfs_recursive(root);
dfs_stack(root);
bfs_queue(root);
//dfs_stack(root);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Commented out code.

for (std::size_t i = 0; i < size; ++i) {
first[i] =
sum(begin(original), end(original), [&](auto const k, auto const el) {
return el * std::exp(i2pi * double(i * k) / double(size));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need to cast to double here?

@Gathros Gathros closed this Jul 17, 2018
@Gathros
Copy link
Contributor Author

Gathros commented Jul 17, 2018

I'm just going to close this pull request and lock it too.

@Gathros Gathros deleted the cppmissingexamples branch July 17, 2018 20:56
@algorithm-archivists algorithm-archivists locked as too heated and limited conversation to collaborators Jul 17, 2018
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Labels
Implementation Edit This provides an edit to an algorithm implementation. (Code and maybe md files are edited.)
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants