PolynomialFunction.cc

Go to the documentation of this file.

#include "casm/basis_set/PolynomialFunction.hh"
 
#include "casm/basis_set/BasisSet.hh"
#include "casm/basis_set/FunctionVisitor.hh"
#include "casm/basis_set/OccupantFunction.hh"
#include "casm/basis_set/Variable.hh"
#include "casm/container/Counter.hh"
#include "casm/container/IsoCounter.hh"
#include "casm/container/MultiCounter.hh"
#include "casm/misc/CASM_Array_math.hh"
#include "casm/symmetry/SymOp.hh"
namespace CASM {
 
PolynomialFunction::PolynomialFunction(
    const std::vector<std::shared_ptr<BasisSet> > &_args)
    : Function(_args), m_coeffs(0) {
  Index i;
  // i counts over subspaces (i.e., BasisSets in init_args)
  for (i = 0; i < _args.size(); i++) {
    // m_subspaces.push_back(Array<Index>(init_args[i].size()));
    // m_sub_sym_reps.push_back(init_args[i].basis_symrep_ID());
    if ((_args[i]->basis_symrep_ID()).empty()) {
      std::cerr << "WARNING: Initializing a PolynomialFunction without knowing "
                   "how to apply symmetry to it. \n"
                << "         Something bad will probably happen; consider this "
                   "your warning!\n";
    }
    if (!_args[i]->size()) {
      std::cerr << "WARNING: In PolynomialFunction constructor, initial "
                   "arguments are ill-defined. Continuing..\n";
      return;
    }
  }
 
  m_coeffs.redefine(num_args());
}
 
//*******************************************************************************************
 
PolynomialFunction::PolynomialFunction(
    const std::vector<std::shared_ptr<BasisSet> > &_args,
    const PolyTrie<double> &_coeffs)
    : Function(_args), m_coeffs(_coeffs) {
  Index i;
  // i counts over subspaces (i.e., BasisSets in init_args)
  for (i = 0; i < _args.size(); i++) {
    // m_subspaces.push_back(Array<Index>(init_args[i].size()));
    // m_sub_sym_reps.push_back(init_args[i].basis_symrep_ID());
    if ((_args[i]->basis_symrep_ID()).empty()) {
      std::cerr << "WARNING: Initializing a PolynomialFunction without knowing "
                   "how to apply symmetry to it. \n"
                << "         Something bad will probably happen; consider this "
                   "your warning!\n";
    }
    if (!_args[i]->size()) {
      std::cerr << "WARNING: In PolynomialFunction constructor, initial "
                   "arguments are ill-defined. Continuing..\n";
      return;
    }
  }
 
  if (m_coeffs.depth() != num_args()) {
    std::cerr << "CRITICAL ERROR: PolynomialFunction initialized with "
                 "coefficients that are incompatible with its argument list.\n"
              << "                Exiting...\n";
    exit(1);
  }
}
 
//*******************************************************************************************
 
PolynomialFunction::PolynomialFunction(const PolynomialFunction &LHS,
                                       const PolynomialFunction &RHS)
    : m_coeffs(0) {
  Index i, j;
  // RHS_ind tracks how exponents of RHS map into resulting product
  Array<Index> RHS_ind;
  for (i = 0; i < LHS.m_argument.size(); i++) {
    m_argument.push_back(LHS.m_argument[i]);
  }
  m_arg2sub = LHS.m_arg2sub;
  m_arg2fun = LHS.m_arg2fun;
  Index linear_ind;
 
  // check each argument set from RHS to see if it already exists in m_argument
  for (j = 0; j < RHS.m_argument.size(); j++) {
    bool add_subspace(true);
    linear_ind = 0;
    for (i = 0; i < LHS.m_argument.size() && add_subspace; i++) {
      if (LHS.m_argument[i]->compare(*RHS.m_argument[j])) {
        add_subspace = false;
        break;
      }
      linear_ind += LHS.m_argument[i]->size();
    }
 
    // RHS.m_argument[j] does not already exist in the produt, so we add it
    if (add_subspace) {
      m_argument.push_back(RHS.m_argument[j]);
      for (i = 0; i < RHS.m_argument[j]->size(); i++) {
        m_arg2sub.push_back(m_argument.size() - 1);
        m_arg2fun.push_back(i);
      }
    }
 
    // linear_ind is the index in (*this) that corresponds to argument
    // (*RHS.m_argument[j])[0]
    RHS_ind.append(Array<Index>::sequence(
        linear_ind, linear_ind + RHS.m_argument[j]->size() - 1));
  }
 
  m_coeffs.redefine(num_args());
  Array<Index> LHS_ind(num_args()), tot_ind(num_args());
 
  PolyTrie<double>::const_iterator LHS_it(LHS.m_coeffs.begin()),
      LHS_end(LHS.m_coeffs.end());
  double t_coeff;
  for (; LHS_it != LHS_end; ++LHS_it) {
    for (i = 0; i < LHS_it.key().size(); i++) {
      LHS_ind[i] = LHS_it.key()[i];
    }
    PolyTrie<double>::const_iterator RHS_it(RHS.m_coeffs.begin()),
        RHS_end(RHS.m_coeffs.end());
    for (; RHS_it != RHS_end; ++RHS_it) {
      t_coeff = (*RHS_it) * (*LHS_it);
      if (almost_zero(t_coeff)) continue;
 
      tot_ind = LHS_ind;
      for (i = 0; i < RHS_it.key().size(); i++) {
        tot_ind[RHS_ind[i]] += RHS_it.key()[i];
      }
      m_coeffs.at(tot_ind) += t_coeff;
    }
  }
  return;
}
 
//*******************************************************************************************
 
PolynomialFunction::PolynomialFunction(const PolynomialFunction &RHS)
    : Function(RHS), m_coeffs(RHS.m_coeffs) {}
 
//*******************************************************************************************
 
PolynomialFunction::PolynomialFunction(const PolynomialFunction &RHS,
                                       const PolyTrie<double> &_coeffs)
    : Function(RHS), m_coeffs(_coeffs) {
  if (m_coeffs.depth() != RHS.m_coeffs.depth()) {
    std::cerr << "WARNING: In PolynomialFunction::PolynomialFunction(const "
                 "PolynomialFunction&, const PolyTrie<double>&),\n"
              << "         the new PolyTrie is incompatible with with the "
                 "number of arguments. Initializing to zero instead.\n";
    m_coeffs.redefine(RHS.m_coeffs.depth());
  }
}
 
//*******************************************************************************************
 
Function *PolynomialFunction::copy() const {
  return new PolynomialFunction(*this);
}
 
//*******************************************************************************************
 
Function *PolynomialFunction::copy(const PolyTrie<double> &new_coeffs) const {
  return new PolynomialFunction(*this, new_coeffs);
}
 
//*******************************************************************************************
 
bool PolynomialFunction::_accept(const FunctionVisitor &visitor,
                                 BasisSet const *home_basis_ptr /*=nullptr*/) {
  return visitor.visit(*this, home_basis_ptr);
}
 
//*******************************************************************************************
 
bool PolynomialFunction::_accept(
    const FunctionVisitor &visitor,
    BasisSet const *home_basis_ptr /*=nullptr*/) const {
  return visitor.visit(*this, home_basis_ptr);
}
 
//*******************************************************************************************
bool PolynomialFunction::depends_on(const Function *test_func) const {
  Index sub_ind(0), arg_ind(0), linear_ind(0);
  for (sub_ind = 0; sub_ind < m_argument.size(); sub_ind++) {
    arg_ind = m_argument[sub_ind]->find(const_cast<Function *const>(test_func));
    linear_ind += arg_ind;
    if (arg_ind < m_argument[sub_ind]->size()) {
      break;
    }
  }
  if (sub_ind == m_argument.size()) return false;
 
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (almost_zero(*it)) continue;
 
    if (it.key()[linear_ind]) return true;
  }
  return false;
}
//*******************************************************************************************
bool PolynomialFunction::is_zero() const {
  // Could check to see if arguments are zero.  For now, assume this is done at
  // time of construction
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (!almost_zero(*it)) return false;
  }
  return true;
}
 
//*******************************************************************************************
void PolynomialFunction::small_to_zero(double tol) {
  m_coeffs.prune_zeros(tol);
}
 
//*******************************************************************************************
Index PolynomialFunction::num_terms() const {
  Index np(0);
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (!almost_zero(*it)) np++;
  }
  return np;
}
 
//*******************************************************************************************
 
double PolynomialFunction::leading_coefficient() const {
  Index t_ind;
  return leading_coefficient(t_ind);
}
 
//*******************************************************************************************
 
double PolynomialFunction::leading_coefficient(Index &index) const {
  index = 0;
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (!almost_zero(*it)) return *it;
    index++;
  }
  return 0.0;
}
 
//*******************************************************************************************
 
double PolynomialFunction::get_coefficient(Index i) const {
  Index index = 0;
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (!almost_zero(*it)) {
      if (index == i) return *it;
 
      index++;
    }
  }
  return 0.0;
}
 
//*******************************************************************************************
 
void PolynomialFunction::make_formula() const {
  PolyTrie<double> tcoeffs(m_coeffs);
  tcoeffs.sort_leaves(ComparePTLeaf::CustomOrder());
  m_formula.clear();
  m_tex_formula.clear();
  // std::cout << "Making PolynomialFunction Formula:\n";
  // tcoeffs.print_sparse(std::cout);
  // std::cout << std::endl;
  std::stringstream tformula, ttex;
  Index np;
  bool is_zero(true);
  Array<Array<Index> > unique_product;
  Array<double> prefactor;
 
  PolyTrie<double>::const_iterator it(tcoeffs.begin()), it_end(tcoeffs.end());
  for (; it != it_end; ++it) {
    if (almost_zero(*it)) continue;
 
    unique_product.push_back(it.key());
    prefactor.push_back(*it);
    is_zero = false;
  }
 
  // Comment out following block to turn off monomial sorting
  // Array<Index> iperm;
  // unique_product.sort(iperm);
  // prefactor.permute(iperm);
  // \end monomial sort
 
  if (is_zero) {
    m_formula = "0";
    m_tex_formula = "0";
    return;
  }
 
  double func_scale(prefactor[0]);
  if (!almost_zero(func_scale - 1)) {
    if (almost_zero(func_scale + 1)) {
      ttex << '-';
    } else {
      ttex << irrational_to_tex_string(func_scale, num_args() * num_args());
    }
  }
 
  if (unique_product.size() > 1) {
    tformula << '(';
    ttex << '(';
  }
 
  for (np = 0; np < unique_product.size(); np++) {
    if (np > 0 && prefactor[np] > 0.0) {
      tformula << '+';
    }
    if (almost_zero(prefactor[np] + 1)) {
      tformula << '-';
    } else if (!almost_zero(prefactor[np] - 1)) {
      tformula << prefactor[np] << " * ";
    }
 
    if (np > 0 && prefactor[np] / func_scale > 0.0) {
      ttex << '+';
    }
    if (almost_zero(prefactor[np] / func_scale + 1)) {
      ttex << '-';
    } else if (!almost_zero(prefactor[np] / func_scale - 1)) {
      ttex << irrational_to_tex_string(prefactor[np] / func_scale,
                                       num_args() * num_args());
    }
 
    // Loop over arguments of unique_product[np]
    int tot_pow(0);
    for (Index linear_ind = 0; linear_ind < unique_product[np].size();
         linear_ind++) {
      if (!unique_product[np][linear_ind]) continue;
      if (tot_pow > 0) {
        tformula << " * ";
      }
      tot_pow += unique_product[np][linear_ind];
 
      if (unique_product[np][linear_ind] > 1) {
        tformula << "pow(" << _argument(linear_ind)->formula() << ", "
                 << unique_product[np][linear_ind] << ")";
        ttex << _argument(linear_ind)->tex_formula() << "^{"
             << unique_product[np][linear_ind] << "} ";
      } else {
        tformula << _argument(linear_ind)->formula();
        ttex << _argument(linear_ind)->tex_formula();
      }
    }
    if (tot_pow == 0 && almost_zero(std::abs(prefactor[np] / func_scale) - 1)) {
      tformula << '1';
      ttex << '1';
    }
  }
 
  if (unique_product.size() > 1) {
    tformula << ')';
    ttex << ')';
  }
  m_tex_formula = ttex.str();
  m_formula = tformula.str();
  // std::cout << "Formula is " << m_formula << '\n';
  // std::cout << "TeX Formula is " << m_tex_formula << '\n';
}
 
//*******************************************************************************************
 
std::set<Index> PolynomialFunction::dof_IDs() const {
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
 
  if (it == it_end) return {};
  std::vector<Index> check(m_coeffs.depth(), false);
  for (; it != it_end; ++it) {
    if (almost_zero(*it)) continue;
 
    for (Index i = 0; i < it.key().size(); i++)
      check[i] = check[i] || it.key()[i];
  }
 
  std::set<Index> result;
  Index last = -1;
  for (Index i = 0; i < check.size(); i++) {
    if (check[i] && last != m_arg2sub[i]) {
      last = m_arg2sub[i];
      std::vector<Index> tres = argument_bases()[last]->dof_IDs();
      result.insert(tres.begin(), tres.end());
    }
  }
  return result;
}
 
//*******************************************************************************************
 
void PolynomialFunction::fill_dispatch_table() {
  Function::inner_prod_table[sclass_ID()][sclass_ID()] =
      new BasicPolyPolyScalarProd();
  Function::operation_table[sclass_ID()][sclass_ID()] = new PolyPolyOperation();
 
  // IMPORTANT: Do
  //       Function::operation_table[OTHER::sclass_ID()][sclass_ID()] = nullptr;
  // Before
  //       Function::operation_table[sclass_ID()][OTHER::sclass_ID()] =
  //       Whatever;
 
  Function::operation_table[Variable::sclass_ID()][sclass_ID()] = nullptr;
  Function::operation_table[OccupantFunction::sclass_ID()][sclass_ID()] =
      nullptr;
  Function::operation_table[sclass_ID()][Variable::sclass_ID()] =
      new PolyVarOperation();
  Function::operation_table[sclass_ID()][OccupantFunction::sclass_ID()] =
      new PolyOccOperation();
 
  // BasicPolyVarProduct does not yet exist
  // Function::inner_prod_table[sclass_ID()][Variable::sclass_ID()]=new
  // BasicPolyVarProduct();
}
 
//*******************************************************************************************
 
SparseTensor<double> const *PolynomialFunction::get_coeffs() const {
  return nullptr;
}
 
//*******************************************************************************************
 
void PolynomialFunction::scale(double scale_factor) {
  m_formula.clear();
  m_tex_formula.clear();
  refresh_ID();
  m_coeffs *= scale_factor;
}
 
//*******************************************************************************************
 
void PolynomialFunction::_set_arguments(
    const ArgumentContainer &new_arg,
    std::vector<Index> const &compatibility_map) {
  if (new_arg.size() < compatibility_map.size() ||
      compatibility_map.size() != argument_bases().size())
    throw std::runtime_error(
        "In PolynomialFunction::_set_arguments called with incompatible "
        "argument size");
  if (new_arg.empty()) return;
 
  std::vector<Index> ind_perm(num_args());
  std::vector<Index> new_offset;
  std::vector<Index> old_offset(1, 0);
 
  Index new_depth = 0;
  for (Index i = 0; i < new_arg.size(); i++) {
    new_offset.push_back(new_depth);
    new_depth += new_arg[i]->size();
  }
 
  Index old_a = 0;
  for (Index i = 0; i < m_argument.size(); i++) {
    Index new_a = new_offset[compatibility_map[i]];
    for (Index j = 0; j < m_argument[i]->size(); j++) {
      ind_perm[old_a++] = new_a++;
    }
  }
 
  PolyTrie<double> new_coeffs(new_depth);
  Array<Index> new_inds(new_depth, 0);
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    for (Index i = 0; i < ind_perm.size(); ++i)
      new_inds[ind_perm[i]] = it.key()[i];
    new_coeffs(new_inds) = *it;
  }
  m_coeffs.swap(new_coeffs);
}
 
//*******************************************************************************************
Function *PolynomialFunction::_apply_sym(const SymOp &op) {
  // std::cout << "Applying symmetry operation to coefficients:\n";
  // m_coeffs.print_sparse(std::cout);
  // m_coeffs.sort_leaves(ComparePTLeaf::CustomOrder());
  m_formula.clear();
  m_tex_formula.clear();
  refresh_ID();
  PolyTrie<double> t_trie(m_coeffs.depth());
  m_coeffs.swap(t_trie);
 
  PolyTrie<double>::const_iterator it(t_trie.begin()), it_end(t_trie.end());
  for (; it != it_end; ++it) transform_monomial_and_add(*it, it.key(), op);
 
  // std::cout << "\n\n And result is:\n";
  // m_coeffs.print_sparse(std::cout);
  return this;
}
 
//*******************************************************************************************
Function *PolynomialFunction::_apply_sym(
    const SymOp &op, const std::vector<bool> &transform_flags) {
  // std::cout << "Applying symmetry operation to coefficients:\n";
  // m_coeffs.print_sparse(std::cout);
  // m_coeffs.sort_leaves(ComparePTLeaf::CustomOrder());
  m_formula.clear();
  m_tex_formula.clear();
  refresh_ID();
  PolyTrie<double> t_trie(m_coeffs.depth());
  m_coeffs.swap(t_trie);
 
  PolyTrie<double>::const_iterator it(t_trie.begin()), it_end(t_trie.end());
  for (; it != it_end; ++it)
    transform_monomial_and_add_new(*it, it.key(), op, transform_flags);
 
  // std::cout << "\n\n And result is:\n";
  // m_coeffs.print_sparse(std::cout);
  return this;
}
 
//*******************************************************************************************
Function *PolynomialFunction::transform_monomial_and_add(
    double prefactor, const Array<Index> &ind, const SymOp &op) {
  assert(ind.size() == m_coeffs.depth() &&
         "\'ind\' Array is not compatible with PolynomialFunction in "
         "PolynomialFunction::transform_monomial_and_add");
 
  Array<Eigen::MatrixXd const *> rep_mats(op.get_matrix_reps(_sub_sym_reps()));
  Array<Array<Array<Index> > > exp_counter(m_argument.size(),
                                           Array<Array<Index> >());
  Array<Array<Array<Index> > > nz_terms(m_argument.size(),
                                        Array<Array<Index> >());
  Array<Array<Array<double> > > nz_coeffs(m_argument.size(),
                                          Array<Array<double> >());
 
  // nz_terms and nz_coeffs keep track of how each individual variable
  // transforms into a linear combination of vectors exp_counter will be used to
  // expand the resulting multinomial expression by counting over the
  // multinomial terms for example, if x^2 transforms to (x-y)^2/2, the end
  // result will be a PolynomialFunction that encodes x^2/2-x*y+y^2/2
  Index linear_offset(0);
  for (Index ns = 0; ns < m_argument.size(); ns++) {
    for (Index na1 = 0; na1 < m_argument[ns]->size(); na1++) {
      nz_terms[ns].push_back(Array<Index>());
      exp_counter[ns].push_back(Array<Index>());
      nz_coeffs[ns].push_back(Array<double>());
 
      if (!ind[linear_offset + na1]) continue;
 
      for (Index na2 = 0; na2 < m_argument[ns]->size(); na2++) {
        if (rep_mats[ns] && !almost_zero((*rep_mats[ns])(na2, na1))) {
          nz_terms[ns][na1].push_back(linear_offset + na2);
          exp_counter[ns][na1].push_back(0);
          nz_coeffs[ns][na1].push_back((*rep_mats[ns])(na2, na1));
        } else if (!rep_mats[ns] &&
                   na1 == na2) {  // assume identity if no symrep exists
          nz_terms[ns][na1].push_back(linear_offset + na2);
          exp_counter[ns][na1].push_back(0);
          nz_coeffs[ns][na1].push_back(1.0);
        }
      }
      exp_counter[ns][na1][0] = ind[linear_offset + na1];
    }
    linear_offset += m_argument[ns]->size();
  }
  /*
    for(int ns=0; ns<exp_counter.size(); ns++){
    for(int na1=0; na1<exp_counter[ns].size(); na1++){
    for(int na2=0; na2<exp_counter[ns][na1].size(); na2++){
    //std::cout << "exp_counter[" << ns << "][" << na1 << "][" << na2 << "] = "
    << exp_counter[ns][na1][na2] << '\n';
    //std::cout << "nz_terms[" << ns << "][" << na1 << "][" << na2 << "] = " <<
    nz_terms[ns][na1][na2] << '\n';
    //std::cout << "nz_coeffs[" << ns << "][" << na1 << "][" << na2 << "] = " <<
    nz_coeffs[ns][na1][na2] << '\n';
    }
    }
    }
  */
  Array<Index> out_ind(ind.size(), 0);
  double out_coeff;
  Index ns, na1, na2;
  bool cflag(true);
  while (cflag) {
    // std::cout << "Still inside while loop and exp_counter is " << exp_counter
    // << '\n';
    out_coeff = prefactor;
    out_ind.resize(ind.size(), 0);
    for (ns = 0; ns < exp_counter.size(); ns++) {
      for (na1 = 0; na1 < exp_counter[ns].size(); na1++) {
        for (na2 = 0; na2 < exp_counter[ns][na1].size(); na2++) {
          out_ind[nz_terms[ns][na1][na2]] += exp_counter[ns][na1][na2];
          out_coeff *= pow(nz_coeffs[ns][na1][na2], exp_counter[ns][na1][na2]);
        }
        out_coeff *= multinomial_coeff(exp_counter[ns][na1]);
      }
    }
    if (out_ind.sum() != ind.sum()) {
      std::cerr << "WARNING: Starting from " << ind
                << " a portion of the result is at " << out_ind << '\n';
    }
    /*
      if(ind != out_ind || !almost_equal(prefactor, out_coeff)) {
      //std::cout << "Monomial term: " << ind << " : " << prefactor << "\n";
      //std::cout << "Transforms to: " << out_ind << " : " << out_coeff << "
      <--"; if(ind != out_ind || !almost_zero(prefactor + out_coeff))
      //std::cout << "******";
      //std::cout << "\n";
      }*/
 
    // this is where the coefficient gets transforemd
    m_coeffs.at(out_ind) += out_coeff;
 
    // Increment exponent counters
    for (ns = 0; ns < exp_counter.size(); ns++) {
      for (na1 = 0; na1 < exp_counter[ns].size(); na1++) {
        if (exp_counter[ns][na1].size() <= 1) continue;
 
        for (na2 = 0; na2 < exp_counter[ns][na1].size() - 1; na2++) {
          if (exp_counter[ns][na1][na2] == 0) continue;
 
          exp_counter[ns][na1][na2]--;
          exp_counter[ns][na1][na2 + 1]++;
          if (na2 > 0) {
            exp_counter[ns][na1][0] = exp_counter[ns][na1][na2];
            exp_counter[ns][na1][na2] = 0;
          }
          break;
        }
        if (na2 == exp_counter[ns][na1].size() - 1) {
          exp_counter[ns][na1][0] = exp_counter[ns][na1][na2];
          exp_counter[ns][na1][na2] = 0;
        } else {
          break;
        }
      }
      if (na1 < exp_counter[ns].size()) break;
    }
    //\done incrementing exponent counters
 
    cflag = ns < exp_counter.size();
  }
 
  // std::cout << "Transformed PolyTrie: \n";
  // m_coeffs.print_sparse(std::cout);
  // std::cout << '\n';
  return this;
}
 
//*******************************************************************************************
// Improved by incorporating MultiCounter and IsoCounter-- not yet tested
 
Function *PolynomialFunction::transform_monomial_and_add_new(
    double prefactor, const Array<Index> &ind, const SymOp &op,
    const std::vector<bool> &transform_flags) {
  assert(ind.size() == m_coeffs.depth() &&
         "\'ind\' Array is not compatible with PolynomialFunction in "
         "PolynomialFunction::transform_monomial_and_add");
  typedef IsoCounter<Array<Index> > TermCounter;
  typedef MultiCounter<TermCounter> SubspaceCounter;
  typedef MultiCounter<SubspaceCounter> ExpCounter;
 
  assert(transform_flags.size() == m_argument.size());
  Array<Eigen::MatrixXd const *> rep_mats(op.get_matrix_reps(_sub_sym_reps()));
  for (Index i = 0; i < transform_flags.size(); i++) {
    if (!transform_flags[i]) rep_mats[i] = nullptr;
  }
  ExpCounter exp_counter;
  Array<Array<Array<Index> > > nz_terms(m_argument.size(),
                                        Array<Array<Index> >());
  Array<Array<Array<double> > > nz_coeffs(m_argument.size(),
                                          Array<Array<double> >());
 
  // nz_terms and nz_coeffs keep track of how each individual variable
  // transforms into a linear combination of vectors exp_counter will be used to
  // expand the resulting multinomial expression by counting over the
  // multinomial terms for example, if x^2 transforms to (x-y)^2/2, the end
  // result will be a genericpolynomialfunction that encodes x^2/2-x*y+y^2/2
 
  // nz_terms and nz_coeffs is basically a sparse depiction of the symrep
  // matrices.  Instead of constructing them here, they should probably be
  // stored in the symreps themselves and accessed via lazy evaluation.
  Index linear_offset(0);
  for (Index ns = 0; ns < m_argument.size(); ns++) {
    exp_counter.push_back(SubspaceCounter());
    for (Index na1 = 0; na1 < m_argument[ns]->size(); na1++) {
      nz_terms[ns].push_back(Array<Index>());
      nz_coeffs[ns].push_back(Array<double>());
 
      int pow_to_distribute(ind[linear_offset + na1]);
 
      if (pow_to_distribute == 0) {
        exp_counter[ns].push_back(TermCounter());
        continue;
      }
 
      for (Index na2 = 0; na2 < m_argument[ns]->size(); na2++) {
        if (rep_mats[ns] && !almost_zero((*rep_mats[ns])(na2, na1))) {
          // nz_flag = true; // indicate that this variable can transform into
          // something else
          nz_terms[ns][na1].push_back(linear_offset + na2);
          nz_coeffs[ns][na1].push_back((*rep_mats[ns])(na2, na1));
        } else if (!rep_mats[ns] &&
                   na1 == na2) {  // assume identity if no symrep exists
          nz_terms[ns][na1].push_back(linear_offset + na2);
          nz_coeffs[ns][na1].push_back(1.0);
        }
      }
      exp_counter[ns].push_back(
          TermCounter(Array<Index>(nz_terms[ns][na1].size(), 0),
                      Array<Index>(nz_terms[ns][na1].size(), pow_to_distribute),
                      1, pow_to_distribute));
    }
    linear_offset += m_argument[ns]->size();
  }
  /*
    for(Index ns=0; ns<exp_counter.size(); ns++){
    for(Index na1=0; na1<exp_counter[ns].size(); na1++){
    for(Index na2=0; na2<exp_counter[ns][na1].size(); na2++){
    //std::cout << "exp_counter[" << ns << "][" << na1 << "][" << na2 << "] = "
    << exp_counter[ns][na1][na2] << '\n';
    //std::cout << "nz_terms[" << ns << "][" << na1 << "][" << na2 << "] = " <<
    nz_terms[ns][na1][na2] << '\n';
    //std::cout << "nz_coeffs[" << ns << "][" << na1 << "][" << na2 << "] = " <<
    nz_coeffs[ns][na1][na2] << '\n';
    }
    }
    }
  */
  Array<Index> out_ind(ind.size(), 0);
  double out_coeff;
  Index ns, na1, na2;
 
  while (exp_counter.valid()) {
    // std::cout << "Still inside while loop and exp_counter is " << exp_counter
    // << '\n';
    out_coeff = prefactor;
    out_ind.resize(ind.size(), 0);
    for (ns = 0; ns < exp_counter.size(); ns++) {
      for (na1 = 0; na1 < exp_counter[ns].size(); na1++) {
        for (na2 = 0; na2 < exp_counter[ns][na1].size(); na2++) {
          out_ind[nz_terms[ns][na1][na2]] += exp_counter[ns][na1][na2];
          out_coeff *= pow(nz_coeffs[ns][na1][na2], exp_counter[ns][na1][na2]);
        }
        out_coeff *= multinomial_coeff(exp_counter[ns][na1]());
      }
    }
    m_coeffs.at(out_ind) += out_coeff;
 
    exp_counter++;
  }
  return this;
}
 
//*******************************************************************************************
 
double PolynomialFunction::remote_eval() const {
  double t_sum(0.0);
  double t_prod;
  Array<double> arg_states(num_args());
 
  for (Index i = 0; i < num_args(); i++)
    arg_states[i] = _argument(i)->remote_eval();
 
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    t_prod = *it;
 
    for (Index i = 0; i < it.key().size(); i++) {
      t_prod *= pow(arg_states[i], it.key()[i]);
    }
 
    t_sum += t_prod;
  }
 
  return t_sum;
}
 
//*******************************************************************************************
 
double PolynomialFunction::remote_deval(const DoF::RemoteHandle &dvar) const {
  double t_sum(0.0);
  double t_prod;
 
  // Array<bool> arg_calc(m_argument.size(), false);
  Array<double> arg_states(num_args());
  Array<double> arg_derivs(num_args());
 
  for (Index i = 0; i < num_args(); i++) {
    arg_states[i] = _argument(i)->remote_eval();
    arg_derivs[i] = _argument(i)->remote_deval(dvar);
  }
 
  Index const *k_begin, *k_end;
 
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (almost_zero(*it)) continue;
    k_begin = it.key().begin();
    k_end = it.key().end();
 
    // calculate monomial derivative
    for (Index const *i = k_begin; i < k_end; i++) {
      if (*i == 0) continue;
      t_prod = (*it) * arg_derivs[i - k_begin];
 
      for (Index const *j = k_begin; j < k_end; j++) {
        for (Index p = 0; p < (*j) - int(i == j); p++)
          t_prod *= arg_states[j - k_begin];
      }
 
      t_sum += t_prod;
    }
  }
 
  return t_sum;
}
 
//*******************************************************************************************
 
double PolynomialFunction::cache_eval() const {
  double t_sum(0.0);
  double t_prod;
 
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    t_prod = *it;
 
    for (Index i = 0; i < it.key().size(); i++) {
      t_prod *= pow(_arg_eval_cache(i), it.key()[i]);
    }
 
    t_sum += t_prod;
  }
 
  return t_sum;
}
 
//*******************************************************************************************
 
double PolynomialFunction::cache_deval(const DoF::RemoteHandle &dvar) const {
  double t_sum(0.0);
  double t_prod;
 
  // Array<bool> arg_calc(m_argument.size(), false);
  // Array<double> arg_states(num_args());
  // Array<double> arg_derivs(num_args());
 
  // for(Index i = 0; i < num_args(); i++) {
  // arg_states[i] = _argument(i)->remote_eval();
  // arg_derivs[i] = _argument(i)->remote_deval(dvar);
  //}
 
  Index const *k_begin, *k_end;
 
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (almost_zero(*it)) continue;
    k_begin = it.key().begin();
    k_end = it.key().end();
 
    // calculate monomial derivative
    for (Index const *i = k_begin; i < k_end; i++) {
      if (*i == 0) continue;
      t_prod = (*it) * _arg_deval_cache(i - k_begin);
 
      for (Index const *j = k_begin; j < k_end; j++) {
        for (Index p = 0; p < (*j) - int(i == j); p++)
          t_prod *= _arg_eval_cache(j - k_begin);
      }
 
      t_sum += t_prod;
    }
  }
 
  return t_sum;
}
 
//*******************************************************************************************
 
double PolynomialFunction::frobenius_scalar_prod(
    const PolynomialFunction &RHS) const {
  // This assumes that all the arguments of *this and RHS are the same and
  // mutually orthogonal This is not generally the case, but it is unclear how
  // to resolve this in the case of the Frobenius product (unlike an inner
  // product based on integration over some domain).
  double prod_result(0.0), tprod(0.0);
  Array<Index> texp;
  // std::cout << "Scalar product of depth " << m_coeffs.depth() << " and depth
  // " << RHS.m_coeffs.depth() << "\n"; std::cout << "               arg_size "
  // << num_args() << " and arg_size " << RHS.num_args() << "\n";
  PolyTrie<double>::const_iterator it(RHS.m_coeffs.begin()),
      it_end(RHS.m_coeffs.end());
  std::vector<std::vector<std::set<Index> > > indep;
  indep.reserve(m_argument.size());
  for (Index i = 0; i < m_argument.size(); i++)
    indep.push_back(m_argument[i]->independent_sub_bases());
 
  for (; it != it_end; ++it) {
    tprod = (*it) * m_coeffs.get(it.key());
    if (almost_zero(tprod, TOL * TOL)) continue;
    Index l = 0;
    for (Index i = 0; i < indep.size(); i++) {
      for (auto const &indep_set : indep[i]) {
        texp.clear();
        for (Index const &func_ind : indep_set) {
          texp.push_back(it.key()[l + func_ind]);
        }
        tprod /= multinomial_coeff(texp);
      }
      l += m_argument[i]->size();
    }
    prod_result += tprod;
  }
  return prod_result;
}
 
//*******************************************************************************************
 
double PolynomialFunction::box_integral_scalar_prod(
    const PolynomialFunction &RHS, double edge_length) const {
  // This assumes that all the arguments of *this and RHS are the same and
  // mutually orthogonal This is not generally the case, and perhaps we should
  // add a check
 
  // returns integral of (*this) x RHS using limits -edge_length to +edge_length
  // (a box of specified edge_length centered at origin) divided by volume of
  // box (so that '1' is always normalized)
 
  double tprod(0), tintegral(0);
  int texp;
  // volume of the box
  double tvol(pow(edge_length, RHS.num_args()));
 
  if (m_coeffs.depth() != RHS.m_coeffs.depth()) {
    std::cerr
        << "WARNING!!! Attempting to get scalar_product between incompatible "
           "PolynomialFunctions. Assuming that they are orthogonal...\n";
    return 0;
  }
 
  PolyTrie<double>::const_iterator RHS_it(RHS.m_coeffs.begin()),
      RHS_end(RHS.m_coeffs.end());
  PolyTrie<double>::const_iterator LHS_it(m_coeffs.begin()),
      LHS_end(m_coeffs.end());
  if (RHS_it == RHS_end || LHS_it == LHS_end) return 0;
 
  for (; RHS_it != RHS_end; ++RHS_it) {
    if (almost_zero(*RHS_it)) continue;
 
    for (LHS_it = m_coeffs.begin(); LHS_it != LHS_end; ++LHS_it) {
      tintegral = (*RHS_it) * (*LHS_it);
      for (Index i = 0; i < LHS_it.key().size(); i++) {
        texp = LHS_it.key()[i] + RHS_it.key()[i];
 
        // Check to see if texp is odd, in which case integral is zero
        if (texp != 2 * (texp / 2)) {
          tintegral = 0;
          break;
        }
 
        tintegral *= 2 * pow(edge_length / 2, texp + 1) / double(texp + 1);
      }
      tprod += tintegral;
    }
  }
  return tprod / tvol;
}
 
//*******************************************************************************************
 
Function *PolynomialFunction::minus_equals(const PolynomialFunction *RHS) {
  m_coeffs -= RHS->m_coeffs;
  return this;
}
 
//*******************************************************************************************
bool PolynomialFunction::compare(const PolynomialFunction *RHS) const {
  // Function::compare() checks that  arguments are equivalent.
 
  // std::cout << " In compare ... " << std::endl;
  // std::cout << " LHS " <<   std::endl;
  // m_coeffs.print_sparse(std::cout);
  // std::cout << " RHS " <<  std::endl;
  //(RHS->m_coeffs).print_sparse(std::cout);
  return m_coeffs.compare(RHS->m_coeffs, TOL);
}
//*******************************************************************************************
bool PolynomialFunction::prune_zeros() { return m_coeffs.prune_zeros(); }
//*******************************************************************************************
Function *PolynomialFunction::plus_equals(const PolynomialFunction *RHS) {
  m_coeffs += RHS->m_coeffs;
  return this;
}
 
//*******************************************************************************************
Function *PolynomialFunction::apply_sym_coeffs(const SymOp &op,
                                               int dependency_layer) {
  std::vector<bool> transform_flags(m_argument.size());
  int tdep;
  for (Index i = 0; i < m_argument.size(); i++) {
    tdep = (m_argument[i]->dependency_layer()) + 1;
    transform_flags[i] = (tdep == dependency_layer);
    if (dependency_layer < tdep) m_argument[i]->apply_sym(op, dependency_layer);
  }
  return _apply_sym(op, transform_flags);
}
 
//*******************************************************************************************
 
Function *PolynomialFunction::poly_quotient(const Variable *RHS) const {
  Index i = 0;
  for (i = 0; i < num_args(); i++) {
    if (_argument(i)->compare(RHS)) break;
  }
  if (i == num_args()) {
    return nullptr;
  }
 
  PolynomialFunction *tpoly(
      new PolynomialFunction(*this, PolyTrie<double>(m_coeffs.depth())));
 
  Array<Index> new_key;
 
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (!it.key()[i]) continue;
 
    new_key = it.key();
    new_key[i]--;
 
    (tpoly->m_coeffs).set(new_key, *it);
  }
 
  return tpoly;
}
 
//*******************************************************************************************
 
Function *PolynomialFunction::poly_remainder(const Variable *RHS) const {
  Index i = 0;
  for (i = 0; i < num_args(); i++) {
    if (_argument(i)->compare(RHS)) break;
  }
  if (i == num_args()) {
    return copy();
  }
 
  PolynomialFunction *tpoly(
      new PolynomialFunction(*this, PolyTrie<double>(m_coeffs.depth())));
 
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (it.key()[i]) continue;
 
    (tpoly->m_coeffs).set(it.key(), *it);
  }
 
  return tpoly;
}
 
//*******************************************************************************************
 
Function *PolynomialFunction::poly_quotient(const OccupantFunction *RHS) const {
  Index i = 0;
  for (i = 0; i < num_args(); i++) {
    if (_argument(i)->compare(RHS)) break;
  }
  if (i == num_args()) {
    return nullptr;
  }
 
  PolynomialFunction *tpoly(
      new PolynomialFunction(*this, PolyTrie<double>(m_coeffs.depth())));
 
  Array<Index> new_key;
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (!it.key()[i]) continue;
 
    new_key = it.key();
    new_key[i]--;
 
    (tpoly->m_coeffs).set(new_key, *it);
  }
 
  return tpoly;
}
 
//*******************************************************************************************
 
Function *PolynomialFunction::poly_remainder(
    const OccupantFunction *RHS) const {
  Index i = 0;
  for (i = 0; i < num_args(); i++) {
    if (_argument(i)->compare(RHS)) break;
  }
  if (i == num_args()) {
    return copy();
  }
 
  PolynomialFunction *tpoly(
      new PolynomialFunction(*this, PolyTrie<double>(m_coeffs.depth())));
 
  PolyTrie<double>::const_iterator it(m_coeffs.begin()), it_end(m_coeffs.end());
  for (; it != it_end; ++it) {
    if (it.key()[i]) continue;
 
    (tpoly->m_coeffs).set(it.key(), *it);
  }
 
  return tpoly;
}
 
//*******************************************************************************************
 
double BasicPolyPolyScalarProd::dot(Function const *LHS,
                                    Function const *RHS) const {
  PolynomialFunction const *PLHS(static_cast<PolynomialFunction const *>(LHS)),
      *PRHS(static_cast<PolynomialFunction const *>(RHS));
 
  return PLHS->frobenius_scalar_prod(*PRHS);
}
 
//*******************************************************************************************
bool PolyPolyOperation::compare(Function const *LHS,
                                Function const *RHS) const {
  PolynomialFunction const *PLHS(static_cast<PolynomialFunction const *>(LHS));
  PolynomialFunction const *PRHS(static_cast<PolynomialFunction const *>(RHS));
  return PLHS->compare(PRHS);
}
 
//*******************************************************************************************
Function *PolyPolyOperation::multiply(Function const *LHS,
                                      Function const *RHS) const {
  PolynomialFunction const *PLHS(static_cast<PolynomialFunction const *>(LHS));
  PolynomialFunction const *PRHS(static_cast<PolynomialFunction const *>(RHS));
 
  return new PolynomialFunction(*PLHS, *PRHS);
}
 
//*******************************************************************************************
 
Function *PolyPolyOperation::multiply_by(Function *LHS,
                                         Function const *RHS) const {
  std::cerr << "WARNING: In-place multiplication of one PolynomialFunction "
               "with another is not yet implemented!!\n"
            << "         Exiting...\n";
  exit(1);
  return nullptr;
}
 
//*******************************************************************************************
Function *PolyPolyOperation::subtract(Function const *LHS,
                                      Function const *RHS) const {
  PolynomialFunction *PLHS_copy(static_cast<PolynomialFunction *>(LHS->copy()));
  PolynomialFunction const *PRHS(static_cast<PolynomialFunction const *>(RHS));
  PLHS_copy->minus_equals(PRHS);
 
  PLHS_copy->refresh_ID();
 
  return PLHS_copy;
}
 
//*******************************************************************************************
Function *PolyPolyOperation::subtract_from(Function *LHS,
                                           Function const *RHS) const {
  PolynomialFunction *PLHS(static_cast<PolynomialFunction *>(LHS));
  PolynomialFunction const *PRHS(static_cast<PolynomialFunction const *>(RHS));
  PLHS->minus_equals(PRHS);
 
  PLHS->refresh_ID();
 
  return PLHS;
}
 
//*******************************************************************************************
Function *PolyPolyOperation::add(Function const *LHS,
                                 Function const *RHS) const {
  PolynomialFunction *PLHS_copy(static_cast<PolynomialFunction *>(LHS->copy()));
  PolynomialFunction const *PRHS(static_cast<PolynomialFunction const *>(RHS));
  PLHS_copy->plus_equals(PRHS);
 
  PLHS_copy->refresh_ID();
 
  return PLHS_copy;
}
 
//*******************************************************************************************
Function *PolyPolyOperation::add_to(Function *LHS, Function const *RHS) const {
  PolynomialFunction *PLHS(static_cast<PolynomialFunction *>(LHS));
  PolynomialFunction const *PRHS(static_cast<PolynomialFunction const *>(RHS));
  PLHS->plus_equals(PRHS);
 
  PLHS->refresh_ID();
 
  return PLHS;
}
 
//*******************************************************************************************
 
Function *PolyVarOperation::poly_quotient(Function const *LHS,
                                          Function const *RHS) const {
  return static_cast<PolynomialFunction const *>(LHS)->poly_quotient(
      static_cast<Variable const *>(RHS));
}
 
//*******************************************************************************************
 
Function *PolyVarOperation::poly_remainder(Function const *LHS,
                                           Function const *RHS) const {
  return static_cast<PolynomialFunction const *>(LHS)->poly_remainder(
      static_cast<Variable const *>(RHS));
}
 
//*******************************************************************************************
 
Function *PolyOccOperation::poly_quotient(Function const *LHS,
                                          Function const *RHS) const {
  return static_cast<PolynomialFunction const *>(LHS)->poly_quotient(
      static_cast<OccupantFunction const *>(RHS));
}
 
//*******************************************************************************************
 
Function *PolyOccOperation::poly_remainder(Function const *LHS,
                                           Function const *RHS) const {
  return static_cast<PolynomialFunction const *>(LHS)->poly_remainder(
      static_cast<OccupantFunction const *>(RHS));
}
 
}  // namespace CASM