C++ expression template example for mixing sparse (indices) and dense vector

functor_template.h

#include <iostream>
#include <vector>
#include <set>
#include <initializer_list>

template<typename E>
struct ExpressionFunctor {
  double operator[](size_t i) const { return static_cast<E const&>(*this)[i]; }
  size_t size()               const { return static_cast<E const&>(*this).size(); }

  operator E& () { return static_cast<E&>(*this); }
  operator const E& () const { return static_cast<const E&>(*this); }
};

struct Vec : ExpressionFunctor<Vec> {
  std::vector<double> elems_;
  double operator[](size_t i) const { return elems_[i]; }
  double &operator[](size_t i)      { return elems_[i]; }
  size_t size() const               { return elems_.size(); }

  Vec(size_t n) : elems_(n) {}
  Vec(std::initializer_list<double> init) {
    for(auto i:init)
      elems_.push_back(i);
  }
};

struct Indices : ExpressionFunctor<Indices> {
  std::set<int> elems_;
  size_t size() const               { return elems_.size(); }

  Indices(std::initializer_list<int> init) {
    for(auto i:init)
      elems_.insert(i);
  }
  // check if has index in indices
  inline bool Has(size_t n) const {
    return elems_.find(n) != elems_.end();
  }
};

template<typename E1, typename E2>
struct SumExp : ExpressionFunctor<SumExp<E1, E2>> {
  E1 const& e1_;
  E2 const& e2_;

  SumExp(E1 const& u, E2 const& v) : e1_(u), e2_(v) {
    std::cout << "construct VecSum" << std::endl;
  }

  // evaluation
  double operator[](size_t i) const { return e1_[i] + e2_[i]; }
  size_t size()               const { return e2_.size(); }
};

template<typename E1, typename IDX, typename E2>
struct SparseAddVec : ExpressionFunctor<SparseAddVec<E1, IDX, E2>> {
  // select using index from v2 to add to v1
  E1  const& e1_;
  IDX const& idx_;
  E2  const& e2_;

  SparseAddVec(E1 const& e1, IDX const& idx, E2 const& e2) :
    e1_(e1), idx_(idx), e2_(e2) {
    std::cout << "construct SparseAddVec" << std::endl;
  }

  // evaluation
  double operator[](size_t i) const {
    if (idx_.Has(i)) {
      return e1_[i] + e2_[i];
    } else {
      return e1_[i]; 
    }
  }
  size_t size()               const { return e1_.size(); }
};

Some code for testing:

#include "functor_template.h"

int main() {
  Vec v1({1,2,3,4});
  Vec v2({1,2,3,4});
  Vec v3({2,2,3,4});
  Indices i({0,2});

  SumExp<Vec, Vec> s1(v1, v2);
  SumExp<SumExp<Vec, Vec>, Vec> s2(s1, v3);
  for (int i = 0; i < s2.size(); ++i)
    std::cout << s2[i] << std::endl;
  
  SparseAddVec<SumExp<Vec, Vec>, Indices, Vec> s3(s1, i, v3);
  for (int i = 0; i < s3.size(); ++i)
    std::cout << s3[i] << std::endl;

  return 0;
}