Commit aeb65415 authored by Nguyen, Thien Minh's avatar Nguyen, Thien Minh
Browse files

Update exp-val calculation based on itensor guidance 



The previous one can cause memory consumption issue. Hence, need to gradually contract the exp-val network.

Signed-off-by: default avatarThien Nguyen <nguyentm@ornl.gov>
parent 55543834

tnqvm/visitors/itensor/mps/ITensorMPSVisitor.cpp

+64 −19
Original line number Diff line number Diff line
@@ -176,28 +176,71 @@ double ITensorMPSVisitor::compute_expectation_z( 
    auto Z_op = singleQubitTensor(getSiteIndex(site_idx), Z(in_measureBits[0]));

    return itensor::eltC(bra * Z_op * ket).real();

  } else {

    auto hamOp = xacc::container::contains(in_measureBits, 0)

                     ? singleQubitTensor(getSiteIndex(1), Z(0))

                     : singleQubitTensor(getSiteIndex(1), Identity(0));



    for (size_t i = 2; i <= m_buffer->size(); ++i) {

      if (xacc::container::contains(in_measureBits, i - 1)) {

        Z obsGate(i - 1);

        auto Op = singleQubitTensor(getSiteIndex(i), obsGate);

        hamOp *= Op;

    // Follow the recipe here:

    // https://www.itensor.org/docs.cgi?vers=cppv3&page=tutorials/correlations

    //'gauge' the MPS to site i

    // any 'position' between i and j, inclusive, would work here

    auto sortedMeasureBits = in_measureBits;

    std::sort(sortedMeasureBits.begin(), sortedMeasureBits.end());

    const auto site_idx = sortedMeasureBits[0] + 1;

    m_mps.position(site_idx);



    auto &psi = m_mps;

    // Create the bra/dual version of the MPS psi

    auto psidag = itensor::dag(m_mps);



    // Prime the link indices to make them distinct from

    // the original ket links

    psidag.prime("Link");

    auto i = site_idx;

    // Handle first qubit: either has measure or not

    auto C = [&]() {

      Z obsGate(site_idx - 1);

      auto op_i = singleQubitTensor(getSiteIndex(site_idx), obsGate);

      if (site_idx != 1) {

        // No measure at q[0]: need to get left link:

        auto li_1 = itensor::leftLinkIndex(psi, site_idx);

        auto Cval = itensor::prime(psi(i), li_1) * op_i;

        Cval *= itensor::prime(psidag(i), "Site");

        return Cval;

      } else {

        Identity obsGate(i - 1);

        auto Op = singleQubitTensor(getSiteIndex(i), obsGate);

        hamOp *= Op;

        auto Cval = psi(i) * op_i;

        Cval *= itensor::prime(psidag(i), "Site");

        return Cval;

      }

    }();

    for (int obs_id = 1; obs_id < sortedMeasureBits.size() - 1; ++obs_id) {

      auto j = sortedMeasureBits[obs_id] + 1;

      for (int k = i + 1; k < j; ++k) {

        C *= psi(k);

        C *= psidag(k);

      }

      Z obsGateJ(sortedMeasureBits[obs_id]);

      auto op_j = singleQubitTensor(getSiteIndex(j), obsGateJ);

      C *= psi(j) * op_j;

      C *= prime(psidag(j), "Site");

      i = j;

    }

    // Last measure bit:

    i = sortedMeasureBits[sortedMeasureBits.size() - 2] + 1;

    auto j = sortedMeasureBits[sortedMeasureBits.size() - 1] + 1;

    for (int k = i + 1; k < j; ++k) {

      C *= psi(k);

      C *= psidag(k);

    }

    Z obsGateJ(j - 1);

    auto op_j = singleQubitTensor(getSiteIndex(j), obsGateJ);

    if (j < m_buffer->size()) {

      auto lj = itensor::rightLinkIndex(psi, j);

      C *= prime(psi(j), lj) * op_j;

      C *= prime(psidag(j), "Site");

    } else {

      C *= psi(j) * op_j;

      C *= prime(psidag(j), "Site");

    }



    auto bond_ket = m_mps(1);

    for (size_t i = 2; i <= m_buffer->size(); ++i) {

      bond_ket *= m_mps(i);

    }

    auto bond_bra = itensor::dag(itensor::prime(bond_ket, "Site"));

    return itensor::eltC(bond_bra * hamOp * bond_ket).real();

    const auto result = itensor::eltC(C);

    // std::cout << "Result = " << result << "\n";

    return result.real();

  }

}
@@ -232,8 +275,10 @@ void ITensorMPSVisitor::finalize() { 
  }



  // Assign the exp-val-z (single Composite mode)

  if (!m_measureBits.empty()) {

    m_buffer->addExtraInfo("exp-val-z", compute_expectation_z(m_measureBits));

  }

}



void ITensorMPSVisitor::applySingleQubitGate(xacc::Instruction &in_gate) {

  assert(in_gate.bits().size() == 1);