VerticalContractor.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef __VERTICAL_CONTRACTOR_H_
#define __VERTICAL_CONTRACTOR_H_ 1
 
#include "BaseContractor.h"
 
#include <boost/container_hash/hash.hpp>
 
namespace TensorNetworks {
 
class VerticalContractor : public BaseContractor {
 public:
  double Contract(const TensorNetwork &network, Types::qubit_t qubit) override {
    std::vector<Eigen::Index> keys;
    std::unordered_map<Eigen::Index, Eigen::Index> keysKeys;
 
    TensorsMap tensors =
        InitializeTensors(network, qubit, keys, keysKeys, false);
 
    const auto &qubitGroup = network.GetQubitGroup(qubit);
    const std::set<Types::qubit_t> qubitGroupSet(qubitGroup.begin(),
                                                 qubitGroup.end());
 
    using TensorPair = std::pair<Eigen::Index, Eigen::Index>;
    std::unordered_set<TensorPair, boost::hash<TensorPair>> visitedPairs;
 
    auto qit = qubitGroupSet.begin();
 
    // while there is more than one tensor...
    while (tensors.size() > 1) {
      if (qit == qubitGroupSet.end())
        throw std::out_of_range("The contraction qubit is out of range.");
 
      const auto q = *qit;
 
      bool found = false;
 
      Eigen::Index tensor1IdBest;
      Eigen::Index tensor2IdBest;
      Eigen::Index resultRankBest;
      Eigen::Index bestCost;
 
      for (auto tensorIt = tensors.begin(); tensorIt != tensors.end();
           ++tensorIt) {
        const auto &tensor = tensorIt->second;
        Eigen::Index tensor1Id = tensor->GetId();
 
        for (Eigen::Index qi = 0;
             qi < static_cast<Eigen::Index>(tensor->qubits.size()); ++qi) {
          if (tensor->qubits[qi] == q) {
            const Eigen::Index tensor2Id = tensor->connections[qi];
 
            if (tensor2Id != TensorNode::NotConnected) {
              auto t1 = tensor1Id;
              auto t2 = tensor2Id;
              if (t1 > t2) std::swap(t1, t2);
              if (visitedPairs.find(std::make_pair(t1, t2)) !=
                  visitedPairs.end())
                continue;
              else
                visitedPairs.insert(std::make_pair(t1, t2));
 
              const Eigen::Index resultRank =
                  GetResultRank(tensors[tensor1Id], tensors[tensor2Id]);
 
              // const Eigen::Index Cost = resultRank -
              // std::max(tensor->GetRank(), tensors[tensor2Id]->GetRank());
              const Eigen::Index Cost =
                  resultRank -
                  (tensor->GetRank() + tensors[tensor2Id]->GetRank()) / 2;
 
              if (!found || resultRank < resultRankBest ||
                  (resultRank == resultRankBest && Cost < bestCost)) {
                resultRankBest = resultRank;
                bestCost = Cost;
                tensor1IdBest = tensor1Id;
                tensor2IdBest = tensor2Id;
                found = true;
              }
            }
          }
        }
      }
 
      visitedPairs.clear();
 
      if (!found) {
        ++qit;
        continue;
      }
 
      ContractNodes(qubit, tensors, tensor1IdBest, tensor2IdBest,
                    resultRankBest);
 
      if (resultRankBest == 0) {
        if (tensors.size() == 1 ||
            tensors[tensor1IdBest]->contractsTheNeededQubit)
          return std::real(tensors[tensor1IdBest]->tensor->atOffset(0));
        // erasing this tensor happens because (not the case anymore, it's
        // avoided) the tensor network might be a disjoint one and a subnetwork
        // is contracted that does not contain the needed qubit
        tensors.erase(tensor1IdBest);
      }
    }
 
    return std::real(tensors.begin()->second->tensor->atOffset(0));
  }
 
  std::shared_ptr<ITensorContractor> Clone() const override {
    auto cloned = std::make_shared<VerticalContractor>();
    cloned->maxTensorRank = maxTensorRank;
    cloned->enableMultithreading = enableMultithreading;
    return cloned;
  }
};
 
}  // namespace TensorNetworks
 
#endif  // __VERTICAL_CONTRACTOR_H_