Tensor.h

Go to the documentation of this file.

#ifndef AMPGEN_TENSOR_H
#define AMPGEN_TENSOR_H
#include <memory.h>
#include <stddef.h>
#include <algorithm>
#include <complex>
#include <initializer_list>
#include <iostream>
#include <string>
#include <vector>
#include <utility>
 
#include "AmpGen/Expression.h"
#include "AmpGen/MsgService.h"
#include "AmpGen/Types.h"
 
#define ADD_DEBUG_TENSOR(X, Y)                                                                                                                                 \
  if(Y != nullptr)                                                                                                                                             \
    for(unsigned i = 0; i < Tensor(X).size(); ++i) Y->emplace_back(std::string(#X) + Tensor::coordinates_to_string(Tensor(X).coords(i)), Tensor(X)[i]);
 
#define ADD_DEBUG_TENSOR_NAMED(X, Y, Z)                                                                                                                        \
  if(Y != nullptr)                                                                                                                                             \
    for(unsigned i = 0; i < X.size(); ++i) Y->emplace_back(Z + Tensor::coordinates_to_string(X.coords(i)), X[i]);
 
namespace AmpGen {
  class ASTResolver;
  class TensorProxy;
 
  class Tensor {
  public:
    class Index {
    private:
      std::shared_ptr<int> m_ptr;
      bool m_isUpper;
 
    public:
      bool operator==(const Tensor::Index &other) const { return m_ptr.get() == other.m_ptr.get(); }
      bool operator!=(const Tensor::Index &other) const { return m_ptr.get() != other.m_ptr.get(); }
      bool isUpper() const { return m_isUpper; };
      explicit Index(bool isUpper = false) : m_ptr(std::make_shared<int>()), m_isUpper(isUpper) {}
      Index(const std::shared_ptr<int> &index, bool isUpper = false) : m_ptr(index), m_isUpper(isUpper) {}
      Index operator-() const { return Index(m_ptr, !m_isUpper); }
      friend std::ostream &operator<<(std::ostream &out, const Index &index);
      uint64_t ptr() const { return uint64_t(m_ptr.get()); }
    };
    struct Dim : public std::vector<unsigned> {
      Dim(unsigned a) : std::vector<unsigned>({a}) {}
      Dim(unsigned a, unsigned b) : std::vector<unsigned>({a, b}) {}
      Dim(unsigned a, unsigned b, unsigned c) : std::vector<unsigned>({a, b, c}) {}
      Dim(unsigned a, unsigned b, unsigned c, unsigned d) : std::vector<unsigned>({a, b, c, d}) {}
    };
 
    Tensor();
    explicit Tensor(const std::vector<Expression> &elements);
    explicit Tensor(const std::vector<unsigned> &dim);
    explicit Tensor(const Tensor::Dim &dim) : Tensor(std::vector<unsigned>(dim)) {};
 
    template <class TYPE> Tensor(const std::initializer_list<TYPE> &elements, const std::vector<unsigned> &dim) : m_dim(dim) {
      setupCoordinates();
      for(auto &x : elements) append(x);
    }
 
    template <class TYPE> Tensor(const std::vector<TYPE> &elements, const std::vector<unsigned> &dim) : m_dim(dim) {
      setupCoordinates();
      for(auto &x : elements) append(x);
    }

    Expression &operator[](const unsigned &i);
    Expression &operator[](const std::vector<unsigned> &co);
    const Expression &operator[](const unsigned &i) const;
    const Expression &operator[](const std::vector<unsigned> &co) const;
 
    Expression get(const unsigned &co);
    Expression get(const unsigned &co) const;
    Expression get(const std::vector<unsigned> &_co) const;

 
    Expression &operator()(const unsigned &a) { return Tensor::operator[]({a}); }
    Expression &operator()(const unsigned &a, const unsigned &b) { return Tensor::operator[]({a, b}); }
    const Expression operator()(const unsigned &a) const { return Tensor::operator[]({a}); }
    const Expression &operator()(const unsigned &a, const unsigned &b) const { return Tensor::operator[]({a, b}); }
 
    TensorProxy operator()(const Tensor::Index &a) const;
    TensorProxy operator()(const Tensor::Index &a, const Tensor::Index &b) const;
    TensorProxy operator()(const Tensor::Index &a, const Tensor::Index &b, const Tensor::Index &c) const;
    TensorProxy operator()(const Tensor::Index &a, const Tensor::Index &b, const Tensor::Index &c, const Tensor::Index &d) const;
    TensorProxy operator()(const std::vector<Tensor::Index> &indices) const;
 
    Tensor operator-() const;
 
    void st(const bool simplify = false);
    bool rankMatches(const Tensor &other);
 
    void imposeSymmetry(unsigned indexA, unsigned indexB);
    void imposeSymmetry(std::vector<unsigned> indices);
 
    Tensor Invert() const;
    std::string to_string(const ASTResolver *resolver = nullptr) const;
 
    int metricSgn(const std::vector<unsigned> &coordinates) const;
    int metricSgn(const unsigned &index) const;
    void append(const Expression &expression);
    void append(const real_t &value);
    void append(const complex_t &value);
    void append(const std::string &value);
    void setupCoordinates();
 
    unsigned nDim() const;
    unsigned rank() const;
    unsigned size() const;
    unsigned index(const std::vector<unsigned> &_co) const;
    unsigned symmetrisedIndex(const std::vector<unsigned> &_co) const;
    unsigned nElements() const;
 
    const std::vector<unsigned> coords(const unsigned &index) const;
    const std::vector<unsigned> &dims() const { return m_dim; }
    const std::string dimString() const;
 
    void print(const bool &eval = false) const;
 
    const std::vector<unsigned> &uniqueElements() const { return m_uniqueElements; }
    void operator+=(const Tensor &rhs);
    void operator-=(const Tensor &rhs);
    Tensor conjugate() const;
 
    static std::vector<unsigned> index_to_coordinates(const unsigned &index, const std::vector<unsigned> &dim);
    static unsigned coordinates_to_index(const std::vector<unsigned> &coords, const std::vector<unsigned> &dim);
    static std::string coordinates_to_string(const std::vector<unsigned> &coordinates);
    template <class... ARGS> static std::vector<unsigned> dim(const ARGS &...args) {
      std::vector<unsigned> rt;
      auto up = std::tuple<ARGS...>(args...);
      for_each(up, [&rt](const unsigned &f) { rt.emplace_back(f); });
      return rt;
    }
 
  private:
    std::vector<unsigned> m_dim;
    std::vector<unsigned> m_symmetrisedCoordinates;
    std::vector<unsigned> m_uniqueElements;
    std::vector<Expression> m_elements;
  };

  class TensorProxy {
  public:
    TensorProxy(const Tensor &tensor, const std::vector<Tensor::Index> &indices);
    std::vector<Tensor::Index> indices() const;
    const Tensor &tensor() const;
    Tensor &tensorMutable();
    TensorProxy reorder(const std::vector<Tensor::Index> &indices);
    operator Expression() const;
    operator Tensor() const;
 
  private:
    Tensor m_tensor;
    std::vector<Tensor::Index> m_indices;
  };
 
  class TensorExpression : public IExpression {
  public:
    TensorExpression(const Tensor &tensor);
    std::string to_string(const ASTResolver *resolver) const override;
    void resolve(ASTResolver &resolver) const override;
    complex_t operator()() const override;
    operator Expression() const;
    Tensor tensor() const { return m_tensor; }
    unsigned size() const { return m_tensor.size(); }
 
  private:
    Tensor m_tensor;
  };
 
  Tensor operator+(const Tensor &, const Tensor &);
  Tensor operator-(const Tensor &, const Tensor &);
  Tensor operator/(const Tensor &, const Expression &);
  Tensor operator*(const Expression &, const Tensor &);
  Tensor operator*(const Tensor &, const Expression &);
 
  Tensor operator/(const Tensor &, const double &);
  Tensor operator*(const double &, const Tensor &);
  Tensor operator*(const Tensor &, const double &);
 
  TensorProxy operator*(const TensorProxy &, const TensorProxy &);
  TensorProxy operator+(const TensorProxy &, const TensorProxy &);
  TensorProxy operator-(const TensorProxy &, const TensorProxy &);
 
  TensorProxy operator/(const TensorProxy &, const Expression &);
  TensorProxy operator*(const Expression &, const TensorProxy &);
  TensorProxy operator*(const TensorProxy &, const Expression &);
 
  TensorProxy operator/(const TensorProxy &, const double &);
  TensorProxy operator*(const double &, const TensorProxy &);
  TensorProxy operator*(const TensorProxy &, const double &);
 
  Tensor Identity(const unsigned &rank = 4);
 
  const Tensor LeviCivita(const unsigned &rank = 4);
  Expression dot(const Tensor &A, const Tensor &B);
 
  std::ostream &operator<<(std::ostream &out, const Tensor::Index &index);
} // namespace AmpGen
 
#endif