Add PolynomialBasisElement class, and an implementation of ChebyshevBasisElement

[hongkai-dai]

As mentioned in #13602 , we want to support richer set of polynomial basis. In this PR, I introduce the base class PolynomialBasis, and the derived class ChebyshevBasis. In the future I will add MonomialBasis and probably other basis like Legendre.

I still need to add more functions to ChebyshevBasis, including Differentiation, operator<<, hash function, etc. But this PR is already large (>600 LOC), so I will add these functions in a future PR.

---

This change is 

[jwnimmer-tri]

I've only just started to look at this PR train now.

The //common:symbolic library is included by most of the code in Drake -- most things want the symbolic::Expression scalar type, even if only indirectly. However, the proportion that will also want these supplemental classes seems much lower. It seems like there is false coupling.

The reason that Expression, Formula, Variable, etc. are all in the same cc_library is that they are intimately related -- once you have Expression, to define comparison operators you need Formula as part of the essential algebra. There is a dependency cycle between these core objects (as evidenced by the forward declarations in their header files).

However, the same cycle reasoning does not apply to ChebyshevPolynomial, the new basis classes, etc. I think they could be a distinct physical dependency layer that sits above the basic types, so that users only pay for them in compilation and linking cost when they actually need them. In short, a new cc_library target (or targets) to separate them out of //common:symbolic.

Perhaps even the basic symbolic_polynomial + symbolic_monomial classes should be split into their own library, distinct from the scalar, for the same reasons.

Also of note, the use of namespace symbolic for files not in a directory named symbolic is contrary to GSG. We have grandfathered in the Expression type & files etc. but we should not really be extending the pattern. We should actually be moving symbolic_expression.h into common/symbolic/expression.h to undo the legacy, eventually.

Also of note, the BsplineBasis lives in //math/.... Perhaps that's a good place for some of these classes as well? The //common folder is pretty full already. Though perhaps a symbolic subdirectory would also solve that.

\CC @soonho-tri to weigh in and help guide this architecture.

[soonho-tri]

We have grandfathered in the Expression type & files etc. but we should not really be extending the pattern. We should actually be moving symbolic_expression.h into common/symbolic/expression.h to undo the legacy, eventually.

I'll open a PR for this task.

We can provide multiple targets under common/symbolic. I think we can provide common which includes symbolic::Expression and its dependencies. @jwnimmer-tri , what do you think?

[jwnimmer-tri]

I think #include "drake/common/symbolic.h" and deps = ["@drake//common:symbolic"] are pervasive. It would be nice to leave those alone without any deprecation warnings, at least for now while we iterate on the newer spellings across several PRs.

As for the cc_library names, before voting I'd have to know the contents of the thing we're naming. It's unclear to me so far exactly which files are in the fundamental subset vs add-ons.

[hongkai-dai]

Perhaps even the basic symbolic_polynomial + symbolic_monomial classes should be split into their own library, distinct from the scalar, for the same reasons.

I very much agree. symbolic polynomial /monomial is mostly used in sum-of-squares optimization in MathematicalProgram, we don't have to include it in symbolic target.

I think symbolic_polynomial.h, symbolic_monomial.h, symbolic_monomial_util.h, symbolic_chebyshev_polynomial.h, symbolic_polynomial_basis.h, symbolic_chebyshev_basis.h can all be put in a new separate drake_cc_library.

Reviewable status: needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits)

[hongkai-dai]

@soonho-tri , may I ask you to review this PR?

Reviewable status: needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits)

[soonho-tri]

+@soonho-tri for platform-review. @hongkai-dai , please find a feature-reviewer for this PR.

Regarding the symbolic reorg, I'll open another issue and cc Hongkai and Jeremy.

Reviewable status: LGTM missing from assignee soonho-tri(platform), needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @soonho-tri)

[hongkai-dai]

To make the purpose of this PR more clear, I explained what I plan to do with the PR trains in the github issue #13803.

Reviewable status: LGTM missing from assignee soonho-tri(platform), needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @soonho-tri)

[hongkai-dai]

+@TobiaMarcucci for feature review please, thanks!

Reviewable status: LGTM missing from assignees TobiaMarcucci,soonho-tri(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @soonho-tri and @TobiaMarcucci)

[soonho-tri]

I'll be back after my dinner...

Reviewed 6 of 10 files at r2.
Reviewable status: 18 unresolved discussions, LGTM missing from assignees TobiaMarcucci,soonho-tri(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai, @soonho-tri, and @TobiaMarcucci)

---

common/symbolic_chebyshev_basis.h, line 10 at r2 (raw file):

#include <map>

#include <Eigen/Core>

nit, not used?

---

common/symbolic_chebyshev_basis.h, line 29 at r2 (raw file):

  explicit ChebyshevBasis(const std::map<Variable, int>& var_to_degree_map);

  ~ChebyshevBasis();

nit, => ~ChebyshevBasis() = default

---

common/symbolic_chebyshev_basis.h, line 32 at r2 (raw file):

  /**
   * Compare two ChebyshevBasis in lexicographic order.

nit, Compare => Compares.

---

common/symbolic_chebyshev_basis.h, line 41 at r2 (raw file):

/**
 * Return the product of two Chebyshev basis.

nit, Return => Returns.

---

common/symbolic_chebyshev_basis.h, line 43 at r2 (raw file):

 * Return the product of two Chebyshev basis.
 * Since Tₘ(x) * Tₙ(x) = 0.5 (Tₘ₊ₙ(x) + Tₘ₋ₙ(x)) if m >= n, the product of
 * Chebyshev basis is the weighted sum of several Chebyshev basises. For example

nit, typo basises => bases?

---

common/symbolic_polynomial_basis.h, line 10 at r2 (raw file):

#include <map>

#include <Eigen/Core>

nit, not used?

---

common/symbolic_polynomial_basis.h, line 21 at r2 (raw file):

 * p(x) = ∑ᵢ cᵢ * ϕᵢ(x), where ϕᵢ(x) is the i'th basis, cᵢ is the coefficient
 * of that basis. The most commonly used basis is monomials. For example
 * in polynomial p(x) = 2x₀²x₁ + 3x₀x1 + 2, x₀²x₁, x₀x1 and 1 are all its basis.

nit, x₀x1 => x₀x₁ (two occurrences in this line).

---

common/symbolic_polynomial_basis.h, line 26 at r2 (raw file):

 * Chebyshev polynomial basis p(x) = 2T₂(x₀)T₁(x₁) + 3T₁(x₁) + 2T₂(x₀),
 * T₂(x₀)T₁(x₁),T₁(x₁), and T₂(x₀) are all its basis. This PolynomialBasis class
 * represents a basis ϕᵢ(x). We can think of a polynomial basis as a mapping f

nit, mapping f => mapping?

---

common/symbolic_polynomial_basis.h, line 31 at r2 (raw file):

 * T₂(x₀)T₁(x₁), it can be thought of as a mapping {x₀ -> 2, x₁ -> 1}.
 *
 * Each of the derived class Derived should implement the following functions

How about,

"Each of the derived class, Derived, should implement the following functions"

?

---

common/symbolic_polynomial_basis.h, line 33 at r2 (raw file):

 * Each of the derived class Derived should implement the following functions
 *
 * std::map<Derived, double> operator*(const Derived& A, const Derived&B)

nit, can you itemize them?

 * - std::map<Derived, double> operator*(const Derived& A, const Derived&B)
 * - std::map<Derived, double> Derived::Differentiate(const Variable& var) const;
 * - bool Derived::operator<(const Derived& other) const;

---

common/symbolic_polynomial_basis.h, line 45 at r2 (raw file):

   * Constructs a polynomial basis given the variable and the degree of that
   * variable.
   * @note the degree should be non-negative.

It's better to use @throw std::logic_error ....

---

common/symbolic_polynomial_basis.h, line 50 at r2 (raw file):

  explicit PolynomialBasis(const std::map<Variable, int>& var_to_degree_map);

  virtual ~PolynomialBasis() {}

nit, slightly better to write virtual ~PolynomialBasis() = default.

---

common/symbolic_polynomial_basis.h, line 90 at r2 (raw file):

  // variable to its degree.
  std::map<Variable, int> var_to_degree_map_;
  int total_degree_;

nit, need {} at the end.

---

common/symbolic_polynomial_basis.cc, line 3 at r2 (raw file):

// NOLINTNEXTLINE(build/include): Its header file is included in symbolic.h.

#include <queue>

nit, not needed?

---

common/symbolic_polynomial_basis.cc, line 4 at r2 (raw file):

#include <queue>

nit, we need

#include <algorithm> // for lexicographical_compare
#include <numeric>  // for accumulate
#include <stdexcept> // for logic_error
#include <typeinfo>  // for typeid
#include <utility> // for pair

---

common/symbolic_polynomial_basis.cc, line 54 at r2 (raw file):

bool PolynomialBasis::operator==(const PolynomialBasis& other) const {
  return typeid(*this) == typeid(other) &&
         this->var_to_degree_map() == other.var_to_degree_map();

This does not work with the following example:

   const DerivedBasisB p5({{y_, 2}});
   const DerivedBasisB p6({{x_, 2}});

as it will form y_ == x_ which is a symbolic formula, not a bool.

You need something like:

  if (typeid(*this) != typeid(other)) {
    return false;
  }
  if (var_to_degree_map().size() != other.var_to_degree_map().size()) {
    return false;
  }
  for (auto it1 = var_to_degree_map().begin(),
            it2 = other.var_to_degree_map().begin();
       it1 != var_to_degree_map().end(); ++it1, ++it2) {
    const Variable& var1{it1->first};
    const Variable& var2{it2->first};
    const int exponent1{it1->second};
    const int exponent2{it2->second};
    if (!var1.equal_to(var2) || exponent1 != exponent2) {
      return false;
    }
  }
  return true;

---

common/test/symbolic_polynomial_basis_test.cc, line 89 at r2 (raw file):

}

TEST_F(SymbolicPolynomialBasisTest, operator_equal_not_equal) {

No underscore in GTEST names, operator_equal_not_equal => OperatorEqualNotEqual?

---

common/test/symbolic_polynomial_basis_test.cc, line 125 at r2 (raw file):

}

TEST_F(SymbolicPolynomialBasisTest, less_than) {

Ditto, less_than => LessThan

[soonho-tri]

Reviewable status: 18 unresolved discussions, LGTM missing from assignees TobiaMarcucci,soonho-tri(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai, @soonho-tri, and @TobiaMarcucci)

---

common/symbolic_polynomial_basis.cc, line 54 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

This does not work with the following example:

   const DerivedBasisB p5({{y_, 2}});
   const DerivedBasisB p6({{x_, 2}});

as it will form y_ == x_ which is a symbolic formula, not a bool.

You need something like:

  if (typeid(*this) != typeid(other)) {
    return false;
  }
  if (var_to_degree_map().size() != other.var_to_degree_map().size()) {
    return false;
  }
  for (auto it1 = var_to_degree_map().begin(),
            it2 = other.var_to_degree_map().begin();
       it1 != var_to_degree_map().end(); ++it1, ++it2) {
    const Variable& var1{it1->first};
    const Variable& var2{it2->first};
    const int exponent1{it1->second};
    const int exponent2{it2->second};
    if (!var1.equal_to(var2) || exponent1 != exponent2) {
      return false;
    }
  }
  return true;

in the snippet, please replace exponent with degree..

[soonho-tri]

First pass completed. PTAL.

Reviewed 4 of 10 files at r2.
Reviewable status: 22 unresolved discussions, LGTM missing from assignees TobiaMarcucci,soonho-tri(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai and @TobiaMarcucci)

---

common/symbolic_chebyshev_basis.cc, line 4 at r2 (raw file):

#include "drake/common/symbolic.h"

missing header

#include <vector>
#include <cmath>

---

common/symbolic_chebyshev_basis.cc, line 45 at r2 (raw file):

  int num_common_variables = 0;
  for (const auto& pair_A : a.var_to_degree_map()) {
    if (b.var_to_degree_map().count(pair_A.first) > 0) {

The complexity of this loop is O(nm) where n = |a.var_to_degree_map()| and m = |b.var_to_degree_map()|. I think we can write a O(n + m) algorithm as a and b are using ordered-maps.

---

common/symbolic_chebyshev_basis.cc, line 53 at r2 (raw file):

  // 2^num_common_variables.
  std::vector<std::map<Variable, int>> chebyshev_basis_all(
      std::pow(2, num_common_variables));

I think the use of std::pow is not optimal here because 1) std::pow return a double not an int and 2) there is a possible overflow / truncation when it implicitly does double => int conversion. I think it's better to create and use a helper function which uses the shift-arithmetic operator (e.g. 1 << n) to compute pow(2, n) and checks the value of n to avoid overflow.

---

common/symbolic_chebyshev_polynomial.h, line 121 at r2 (raw file):

/**
 * Evaluate a Chebyshev polynomial at a given value.

nit, Evaluate => Evaluates

[TobiaMarcucci]

Reviewed 10 of 10 files at r2.
Reviewable status: 24 unresolved discussions, LGTM missing from assignees TobiaMarcucci,soonho-tri(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai)

a discussion (no related file):
I've nothing to point out on the coding side, all the functions seem implemented correctly to me. I like the fact that you neglect zero-degree entries in the PolynomialBasis map, it makes life much easier. And even if the implementation of the Chebyshev product ends up being pretty long, I think there is no easy way around it.

I have an observation regarding the wording, similar to the one I wrote in #13803. I think that it isn't fully correct to call phi_i(x) a "basis". My understanding is that the basis is the set {phi_1, ..., phi_n}, whereas each of the functions phi_i(x) is a basis vector. (I think this is also the common way to call them: https://en.wikipedia.org/wiki/Basis_(linear_algebra).) I think that a name like PolynomialBasisVector would be mathematically more rigorous than PolynomialBasis. The same for ChebyshevBasis (=> ChebyshevBasisVector). I don't know however if things end up being too verbose, so I don't necessarily suggest this change. In any case, I'd suggest to be clearer in the docs (see the note I left in symbolic_polynomial_basis.h).

---

---

common/symbolic_polynomial_basis.h, line 18 at r2 (raw file):

namespace symbolic {
/**
 * Each polynomial p(x) can be written as a linear combination of its basis

As I wrote above, I'd suggest to call the functions ϕᵢ basis vectors (or elements of the basis) instead of basis, and to leave the word basis for a collection of basis vectors. So I'd suggest: "a linear combination of its basis" => "a linear combination of the elements of its basis"; "where ϕᵢ(x) is the i'th basis" => "where ϕᵢ(x) is the i'th basis vector"; etc...

[soonho-tri]

Reviewable status: 25 unresolved discussions, LGTM missing from assignees TobiaMarcucci,soonho-tri(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai)

---

common/symbolic_polynomial_basis.cc, line 52 at r2 (raw file):

}

bool PolynomialBasis::operator==(const PolynomialBasis& other) const {

This implementation assumes that all of its derived classes will not have data members other than var_to_degree_map. Can you check if this is the case?

If not, it is better to have a virtual member function bool T::EqualTo(const T& other) and call it here:

    return typeid(*this) == typeid(other) && EqualTo(other);

[hongkai-dai]

Reviewable status: 11 unresolved discussions, LGTM missing from assignees TobiaMarcucci,soonho-tri(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai, @soonho-tri, and @TobiaMarcucci)

a discussion (no related file):

Previously, TobiaMarcucci (Tobia Marcucci) wrote…

I've nothing to point out on the coding side, all the functions seem implemented correctly to me. I like the fact that you neglect zero-degree entries in the PolynomialBasis map, it makes life much easier. And even if the implementation of the Chebyshev product ends up being pretty long, I think there is no easy way around it.

I have an observation regarding the wording, similar to the one I wrote in #13803. I think that it isn't fully correct to call phi_i(x) a "basis". My understanding is that the basis is the set {phi_1, ..., phi_n}, whereas each of the functions phi_i(x) is a basis vector. (I think this is also the common way to call them: https://en.wikipedia.org/wiki/Basis_(linear_algebra).) I think that a name like PolynomialBasisVector would be mathematically more rigorous than PolynomialBasis. The same for ChebyshevBasis (=> ChebyshevBasisVector). I don't know however if things end up being too verbose, so I don't necessarily suggest this change. In any case, I'd suggest to be clearer in the docs (see the note I left in symbolic_polynomial_basis.h).

I agree that PolynomialBasis is not the right word. As mentioned in the github issue comment #13803 (comment), what do you think if we change it to PolynomialBasisElement. Since the basis is a set, and what we implement is an element of this set.

---

---

common/symbolic_chebyshev_basis.cc, line 4 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

missing header

#include <vector>
#include <cmath>

Done.

---

common/symbolic_chebyshev_basis.cc, line 45 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

The complexity of this loop is O(nm) where n = |a.var_to_degree_map()| and m = |b.var_to_degree_map()|. I think we can write a O(n + m) algorithm as a and b are using ordered-maps.

Done. Good call

---

common/symbolic_chebyshev_basis.cc, line 53 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

I think the use of std::pow is not optimal here because 1) std::pow return a double not an int and 2) there is a possible overflow / truncation when it implicitly does double => int conversion. I think it's better to create and use a helper function which uses the shift-arithmetic operator (e.g. 1 << n) to compute pow(2, n) and checks the value of n to avoid overflow.

Done.

---

common/symbolic_polynomial_basis.h, line 31 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

How about,

"Each of the derived class, Derived, should implement the following functions"

?

Done.

---

common/symbolic_polynomial_basis.h, line 45 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

It's better to use @throw std::logic_error ....

Done.

---

common/symbolic_polynomial_basis.cc, line 4 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

nit, we need

#include <algorithm> // for lexicographical_compare
#include <numeric>  // for accumulate
#include <stdexcept> // for logic_error
#include <typeinfo>  // for typeid
#include <utility> // for pair

Thanks! Now I am wondering why the lint test doesn't complain about missing headers. It used to remind me to include these headers explicitly.

---

common/symbolic_polynomial_basis.cc, line 52 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

This implementation assumes that all of its derived classes will not have data members other than var_to_degree_map. Can you check if this is the case?

If not, it is better to have a virtual member function bool T::EqualTo(const T& other) and call it here:

    return typeid(*this) == typeid(other) && EqualTo(other);

Done.

---

common/symbolic_polynomial_basis.cc, line 54 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

in the snippet, please replace exponent with degree..

Good catch! Thanks

---

common/test/symbolic_polynomial_basis_test.cc, line 89 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

No underscore in GTEST names, operator_equal_not_equal => OperatorEqualNotEqual?

Done.

---

common/test/symbolic_polynomial_basis_test.cc, line 125 at r2 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

Ditto, less_than => LessThan

Done.

[TobiaMarcucci]

As said in #13803, I think that both PolynomialBasisElement and PolynomialBasisFunction are valid choices.

[soonho-tri]

Reviewed 6 of 6 files at r3.
Reviewable status: 3 unresolved discussions, LGTM missing from assignees TobiaMarcucci,soonho-tri(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai and @TobiaMarcucci)

---

common/symbolic_chebyshev_basis.cc, line 36 at r3 (raw file):

  } else {
    throw std::invalid_argument("power of 2 overflow");
  }

power_of_2(35) still returns 8 which is not the correct answer. Please consider the following alternative.

int power_of_two(int n) {
  if (n < 0) {
    // handle underflow
  }
  if (n > (sizeof(int) * CHAR_BIT - 2)) {  // need #include <cmath> for CHAR_BIT
    // handle overflow
  }
  return 1 << n;
}

[hongkai-dai]

Reviewable status: 3 unresolved discussions, LGTM missing from assignees TobiaMarcucci,soonho-tri(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai, @soonho-tri, and @TobiaMarcucci)

---

common/symbolic_chebyshev_basis.cc, line 36 at r3 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

power_of_2(35) still returns 8 which is not the correct answer. Please consider the following alternative.

int power_of_two(int n) {
  if (n < 0) {
    // handle underflow
  }
  if (n > (sizeof(int) * CHAR_BIT - 2)) {  // need #include <cmath> for CHAR_BIT
    // handle overflow
  }
  return 1 << n;
}

Done. Thanks!

[soonho-tri]

Reviewed 1 of 1 files at r4.
Reviewable status: 2 unresolved discussions, LGTM missing from assignee TobiaMarcucci, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai and @TobiaMarcucci)

[hongkai-dai]

Reviewable status: 2 unresolved discussions, LGTM missing from assignee TobiaMarcucci, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @soonho-tri and @TobiaMarcucci)

a discussion (no related file):

Previously, hongkai-dai (Hongkai Dai) wrote…

I agree that PolynomialBasis is not the right word. As mentioned in the github issue comment #13803 (comment), what do you think if we change it to PolynomialBasisElement. Since the basis is a set, and what we implement is an element of this set.

Done. I changed the name to PolynomialBasisElement and ChebyshevBasisElement.

---

---

common/symbolic_polynomial_basis.h, line 18 at r2 (raw file):

Previously, TobiaMarcucci (Tobia Marcucci) wrote…

As I wrote above, I'd suggest to call the functions ϕᵢ basis vectors (or elements of the basis) instead of basis, and to leave the word basis for a collection of basis vectors. So I'd suggest: "a linear combination of its basis" => "a linear combination of the elements of its basis"; "where ϕᵢ(x) is the i'th basis" => "where ϕᵢ(x) is the i'th basis vector"; etc...

Done.

[hongkai-dai]

Reviewable status: 3 unresolved discussions, LGTM missing from assignee TobiaMarcucci, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @soonho-tri and @TobiaMarcucci)

---

common/symbolic_chebyshev_basis_element.cc, line 47 at r5 (raw file):

std::map<ChebyshevBasisElement, double> operator*(
    const ChebyshevBasisElement& a, const ChebyshevBasisElement& b) {
  // If variable x shows up in both ChebyshevBasisElement a and b, we know that

@soonho-tri this implementation is long and error-prone. What I really want to do is a cartesian product between some vectors. Say I have

a = T1(x) T2(y)
b = T3(x) T3(y)

I know T1(x) * T3(x) = 0.5 * (T2(x) + T4(x)), T2(y) * T3(y) = 0.5 * (T1(y) * T5(y)). Then to compute their product a * b, I want to have a Cartesian product of
(T2(x), T4(x)) ⊗ (T1(y), T5(y)), which should give me 4 tuples
(T2(x), T1(y))
(T2(x), T5(y))
(T4(x), T1(y))
(T4(x), T5(y))

In python this is simple, we just need to use itertools.product https://docs.python.org/2/library/itertools.html#itertools.product. Unfortunately in C++ this becomes much harder. I cannot use the nested for loop approach either, since the number of for loops is undetermined (depends on how many common variables are shared in a and b.)

Do you know a C++ method that can do cartesian product between several number of vectors? (The number of vectors is unknown until runtime).

[soonho-tri]

Reviewable status: 3 unresolved discussions, LGTM missing from assignee TobiaMarcucci, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai, @soonho-tri, and @TobiaMarcucci)

---

common/symbolic_chebyshev_basis_element.cc, line 47 at r5 (raw file):

Previously, hongkai-dai (Hongkai Dai) wrote…

@soonho-tri this implementation is long and error-prone. What I really want to do is a cartesian product between some vectors. Say I have

a = T1(x) T2(y)
b = T3(x) T3(y)

I know T1(x) * T3(x) = 0.5 * (T2(x) + T4(x)), T2(y) * T3(y) = 0.5 * (T1(y) * T5(y)). Then to compute their product a * b, I want to have a Cartesian product of
(T2(x), T4(x)) ⊗ (T1(y), T5(y)), which should give me 4 tuples
(T2(x), T1(y))
(T2(x), T5(y))
(T4(x), T1(y))
(T4(x), T5(y))

In python this is simple, we just need to use itertools.product https://docs.python.org/2/library/itertools.html#itertools.product. Unfortunately in C++ this becomes much harder. I cannot use the nested for loop approach either, since the number of for loops is undetermined (depends on how many common variables are shared in a and b.)

Do you know a C++ method that can do cartesian product between several number of vectors? (The number of vectors is unknown until runtime).

I see. Mr. Google points to this.

[hongkai-dai]

Reviewable status: 3 unresolved discussions, LGTM missing from assignee TobiaMarcucci, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @soonho-tri and @TobiaMarcucci)

---

common/symbolic_chebyshev_basis_element.cc, line 47 at r5 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

I see. Mr. Google points to this.

Thanks! Hopefully C++ would become more like python with these itertools in the standard library.

[TobiaMarcucci]

Reviewable status: 2 unresolved discussions, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @soonho-tri and @TobiaMarcucci)

[TobiaMarcucci]

Reviewed 1 of 6 files at r3, 8 of 9 files at r5.
Reviewable status: 1 unresolved discussion, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @soonho-tri)

[soonho-tri]

Reviewed 8 of 9 files at r5.
Reviewable status: 1 unresolved discussion, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @hongkai-dai)

---

common/symbolic_chebyshev_basis_element.cc, line 47 at r5 (raw file):

Previously, hongkai-dai (Hongkai Dai) wrote…

Thanks! Hopefully C++ would become more like python with these itertools in the standard library.

Related articles:

• http://conradomiranda.com/projects/libraries/cpp-product-iterator/
• https://github.com/mirandaconrado/product-iterator

If you want, please add TODO and consider opening another PR later.

[hongkai-dai]

Reviewable status: commits need curation (https://drake.mit.edu/reviewable.html#curated-commits)

---

common/symbolic_chebyshev_basis_element.cc, line 47 at r5 (raw file):

Previously, soonho-tri (Soonho Kong) wrote…

Related articles:

• http://conradomiranda.com/projects/libraries/cpp-product-iterator/
• https://github.com/mirandaconrado/product-iterator

If you want, please add TODO and consider opening another PR later.

Thanks Soonho, these articles are helpful.

I think I will stick with the current implementation for now, without adding more dependencies (like product_iterator) to Drake.