Commit 89528309 authored by Mccaskey, Alex's avatar Mccaskey, Alex
Browse files

adding qjit.as_unitary_matrix() to python api 



Signed-off-by: Mccaskey, Alex's avatarAlex McCaskey <mccaskeyaj@ornl.gov>
parent 60cdbb84

python/examples/as_unitary_matrix.py

0 → 100644
+25 −0
Original line number Diff line number Diff line 
from qcor import *

import numpy as np



@qjit

def ansatz(q : qreg, x : float):

    X(q[0])

    Ry(q[1], x)

    CX(q[1], q[0])

        

# Can convert qjit kernels to a unitary matrix form

# I know the ground state is at x = .59

u_mat = ansatz.as_unitary_matrix(qalloc(2), .59)



# Create a Hamiltonian, map it to a numpy matrix

H = -2.1433 * X(0) * X(1) - 2.1433 * \

    Y(0) * Y(1) + .21829 * Z(0) - 6.125 * Z(1) + 5.907

Hmat = H.to_numpy()



# Compute |psi> = U |0>

zero_state = np.array([1., 0., 0., 0.])

final_state = np.dot(u_mat, np.transpose(zero_state))



# Compute E = <psi| H |psi>

energy = np.dot(final_state, np.dot(Hmat,final_state))

print('Energy: ', energy)
 
 No newline at end of file

python/py-qcor.cpp

+23 −1
Original line number Diff line number Diff line
@@ -144,7 +144,9 @@ std::shared_ptr<qcor::Observable> convertToQCOROperator( 


    } else {

      if (keep_fermion) {

        xacc::error("Error - you asked for a qcor::FermionOperator, but this is an OpenFermion QubitOperator.");

        xacc::error(

            "Error - you asked for a qcor::FermionOperator, but this is an "

            "OpenFermion QubitOperator.");

      }

      // this is a qubit  operator

      auto terms = op.attr("terms");
@@ -496,6 +498,7 @@ PYBIND11_MODULE(_pyqcor, m) { 
            qjit.invoke_with_hetmap(name, m);

          },

          "")

          

      .def("extract_composite",

           [](qcor::QJIT &qjit, const std::string name, KernelArgDict args) {

             xacc::HeterogeneousMap m;
@@ -504,6 +507,25 @@ PYBIND11_MODULE(_pyqcor, m) { 
               mpark::visit(vis, item.second);

             }

             return qjit.extract_composite_with_hetmap(name, m);

           })

      .def("internal_as_unitary",

           [](qcor::QJIT &qjit, const std::string name, KernelArgDict args) {

             xacc::HeterogeneousMap m;

             for (auto &item : args) {

               KernelArgDictToHeterogeneousMap vis(m, item.first);

               mpark::visit(vis, item.second);

             }

             auto composite = qjit.extract_composite_with_hetmap(name, m);

             auto n_qubits = composite->nLogicalBits();

             qcor::KernelToUnitaryVisitor visitor(n_qubits);

             InstructionIterator iter(composite);

             while (iter.hasNext()) {

               auto inst = iter.next();

               if (!inst->isComposite() && inst->isEnabled()) {

                 inst->accept(&visitor);

               }

             }

             return visitor.getMat();

           });



  py::class_<qcor::ObjectiveFunction, std::shared_ptr<qcor::ObjectiveFunction>>(

python/qcor.in.py

+7 −1
Original line number Diff line number Diff line
@@ -510,6 +510,12 @@ class qjit(object): 
        """

        return self.extract_composite(*args).nInstructions()

    

    def as_unitary_matrix(self, *args):

        args_dict = {}

        for i, arg_name in enumerate(self.arg_names):

            args_dict[arg_name] = list(args)[i]

        return self._qjit.internal_as_unitary(self.function.__name__, args_dict)

    

    def ctrl(self, *args):

        print('This is an internal API call and will be translated to C++ via the QJIT.\nIt can only be called from within another quantum kernel.')

        exit(1)

python/tests/test_kernel_jit.py

+27 −0
Original line number Diff line number Diff line
@@ -6,6 +6,33 @@ import math as myMathMod 
# Some global variables for testing

MY_PI = 3.1416

class TestKernelJIT(unittest.TestCase):

    def test_as_unitary(self):

        try:

           import numpy as np

        except:

            print('No numpy, cant run test_as_unitary')

            return

        

        @qjit

        def dansatz(q : qreg, x : float):

            X(q[0])

            Ry(q[1], x)

            CX(q[1], q[0])

        

        u_mat = dansatz.as_unitary_matrix(qalloc(2), .59)



        H = -2.1433 * X(0) * X(1) - 2.1433 * \

            Y(0) * Y(1) + .21829 * Z(0) - 6.125 * Z(1) + 5.907

        Hmat = H.to_numpy()



        # Compute |psi> = U |0>

        zero_state = np.array([1., 0., 0., 0.])

        final_state = np.dot(u_mat, np.transpose(zero_state))

        # Compute E = <psi| H |psi>

        energy = np.dot(final_state, np.dot(Hmat,final_state))

        print(energy)

        self.assertAlmostEqual(energy, -1.74, places=1)



    def test_rewrite_decompose(self):

        set_qpu('qpp', {'shots':1024})

        @qjit

runtime/utils/qcor_utils.cpp

+254 −12
Original line number Diff line number Diff line 
#include "qcor_utils.hpp"



#include "qrt.hpp"

#include "xacc.hpp"

#include "xacc_service.hpp"
@@ -32,8 +33,8 @@ std::shared_ptr<qcor::CompositeInstruction> create_ctrl_u() { 
  return std::dynamic_pointer_cast<xacc::CompositeInstruction>(

      xacc::getService<xacc::Instruction>("C-U"));

}

std::shared_ptr<qcor::IRTransformation>

get_transformation(const std::string &transform_type) {

std::shared_ptr<qcor::IRTransformation> get_transformation(

    const std::string &transform_type) {

  if (!xacc::isInitialized())

    xacc::internal_compiler::compiler_InitializeXACC();

  return xacc::getService<xacc::IRTransformation>(transform_type);
@@ -43,8 +44,8 @@ std::shared_ptr<qcor::IRProvider> get_provider() { 
  return xacc::getIRProvider("quantum");

}



std::shared_ptr<qcor::CompositeInstruction>

decompose_unitary(const std::string algorithm, UnitaryMatrix &mat,

std::shared_ptr<qcor::CompositeInstruction> decompose_unitary(

    const std::string algorithm, UnitaryMatrix &mat,

    const std::string buffer_name) {

  std::string local_algorithm = algorithm;

  if (local_algorithm == "z_y_z") local_algorithm = "z-y-z";
@@ -72,10 +73,9 @@ decompose_unitary(const std::string algorithm, UnitaryMatrix &mat, 
  return decomposed;

}



std::shared_ptr<qcor::CompositeInstruction>

decompose_unitary(const std::string algorithm, UnitaryMatrix &mat,

                  const std::string buffer_name,

                  std::shared_ptr<xacc::Optimizer> optimizer) {

std::shared_ptr<qcor::CompositeInstruction> decompose_unitary(

    const std::string algorithm, UnitaryMatrix &mat,

    const std::string buffer_name, std::shared_ptr<xacc::Optimizer> optimizer) {

  auto tmp = xacc::getService<xacc::Instruction>(algorithm);

  auto decomposed = std::dynamic_pointer_cast<CompositeInstruction>(tmp);
@@ -99,4 +99,246 @@ decompose_unitary(const std::string algorithm, UnitaryMatrix &mat, 
}



}  // namespace __internal__

using namespace xacc::quantum;



MatrixXcd X_Mat{MatrixXcd::Zero(2, 2)};

MatrixXcd Y_Mat{MatrixXcd::Zero(2, 2)};

MatrixXcd Z_Mat{MatrixXcd::Zero(2, 2)};

MatrixXcd Id_Mat{MatrixXcd::Identity(2, 2)};

MatrixXcd H_Mat{MatrixXcd::Zero(2, 2)};

MatrixXcd T_Mat{MatrixXcd::Zero(2, 2)};

MatrixXcd Tdg_Mat{MatrixXcd::Zero(2, 2)};

MatrixXcd S_Mat{MatrixXcd::Zero(2, 2)};

MatrixXcd Sdg_Mat{MatrixXcd::Zero(2, 2)};

MatrixXcd CX_Mat{MatrixXcd::Zero(4, 4)};

MatrixXcd CY_Mat{MatrixXcd::Zero(4, 4)};

MatrixXcd CZ_Mat{MatrixXcd::Zero(4, 4)};

MatrixXcd CH_Mat{MatrixXcd::Zero(4, 4)};

MatrixXcd Swap_Mat{MatrixXcd::Zero(4, 4)};



MatrixXcd kroneckerProduct(MatrixXcd &lhs, MatrixXcd &rhs) {

  MatrixXcd result(lhs.rows() * rhs.rows(), lhs.cols() * rhs.cols());

  for (int i = 0; i < lhs.rows(); ++i) {

    for (int j = 0; j < lhs.cols(); ++j) {

      result.block(i * rhs.rows(), j * rhs.cols(), rhs.rows(), rhs.cols()) =

          lhs(i, j) * rhs;

    }

  }

  return result;

}



KernelToUnitaryVisitor::KernelToUnitaryVisitor(size_t in_nbQubits)

    : m_nbQubit(in_nbQubits) {

  static bool static_gate_init = false;

  m_circuitMat = MatrixXcd::Identity(1ULL << in_nbQubits, 1ULL << in_nbQubits);

  if (!static_gate_init) {

    X_Mat << 0.0 + 0.0_i, 1.0 + 0.0_i, 1.0 + 0.0_i, 0.0 + 0.0_i;

    Y_Mat << 0.0 + 0.0_i, -I, I, 0.0 + 0.0_i;

    Z_Mat << 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, -1.0 + 0.0_i;

    H_Mat << 1.0 / std::sqrt(2.0) + 0.0_i, 1.0 / std::sqrt(2.0) + 0.0_i,

        1.0 / std::sqrt(2.0) + 0.0_i, -1.0 / std::sqrt(2.0) + 0.0_i;

    T_Mat << 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, std::exp(I * M_PI / 4.0);

    Tdg_Mat << 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, std::exp(-I * M_PI / 4.0);

    S_Mat << 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, I;

    Sdg_Mat << 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, -I;

    CX_Mat << 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        0.0 + 0.0_i, 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 1.0 + 0.0_i,

        0.0 + 0.0_i;

    CY_Mat << 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        0.0 + 0.0_i, -I, 0.0 + 0.0_i, 0.0 + 0.0_i, I, 0.0 + 0.0_i;

    CZ_Mat << 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        -1.0 + 0.0_i;

    CH_Mat << 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        1.0 / std::sqrt(2.0) + 0.0_i, 1.0 / std::sqrt(2.0) + 0.0_i, 0.0 + 0.0_i,

        0.0 + 0.0_i, 1.0 / std::sqrt(2.0) + 0.0_i,

        -1.0 / std::sqrt(2.0) + 0.0_i;

    Swap_Mat << 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        0.0 + 0.0_i, 1.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 1.0 + 0.0_i,

        0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i, 0.0 + 0.0_i,

        1.0 + 0.0_i;

    static_gate_init = true;

  }

}

void KernelToUnitaryVisitor::visit(Hadamard &h) {

  m_circuitMat = singleQubitGateExpand(H_Mat, h.bits()[0]) * m_circuitMat;

}



void KernelToUnitaryVisitor::visit(CNOT &cnot) {

  MatrixXcd fullMat =

      twoQubitGateExpand(CX_Mat, cnot.bits()[0], cnot.bits()[1]);

  m_circuitMat = fullMat * m_circuitMat;

}



void KernelToUnitaryVisitor::visit(Rz &rz) {

  double theta = InstructionParameterToDouble(rz.getParameter(0));

  MatrixXcd fullMat =

      singleParametricQubitGateExpand(theta, Rot::Z, rz.bits()[0]);

  m_circuitMat = fullMat * m_circuitMat;

}



void KernelToUnitaryVisitor::visit(Ry &ry) {

  double theta = InstructionParameterToDouble(ry.getParameter(0));

  MatrixXcd fullMat =

      singleParametricQubitGateExpand(theta, Rot::Y, ry.bits()[0]);

  m_circuitMat = fullMat * m_circuitMat;

}



void KernelToUnitaryVisitor::visit(Rx &rx) {

  double theta = InstructionParameterToDouble(rx.getParameter(0));

  MatrixXcd fullMat =

      singleParametricQubitGateExpand(theta, Rot::X, rx.bits()[0]);

  m_circuitMat = fullMat * m_circuitMat;

}



void KernelToUnitaryVisitor::visit(X &x) {

  m_circuitMat = singleQubitGateExpand(X_Mat, x.bits()[0]) * m_circuitMat;

}



void KernelToUnitaryVisitor::visit(Y &y) {

  m_circuitMat = singleQubitGateExpand(Y_Mat, y.bits()[0]) * m_circuitMat;

}

void KernelToUnitaryVisitor::visit(Z &z) {

  m_circuitMat = singleQubitGateExpand(Z_Mat, z.bits()[0]) * m_circuitMat;

}

void KernelToUnitaryVisitor::visit(CY &cy) {

  MatrixXcd fullMat = twoQubitGateExpand(CY_Mat, cy.bits()[0], cy.bits()[1]);

  m_circuitMat = fullMat * m_circuitMat;

}

void KernelToUnitaryVisitor::visit(CZ &cz) {

  MatrixXcd fullMat = twoQubitGateExpand(CZ_Mat, cz.bits()[0], cz.bits()[1]);

  m_circuitMat = fullMat * m_circuitMat;

}

void KernelToUnitaryVisitor::visit(Swap &s) {

  MatrixXcd fullMat = twoQubitGateExpand(Swap_Mat, s.bits()[0], s.bits()[1]);

  m_circuitMat = fullMat * m_circuitMat;

}

void KernelToUnitaryVisitor::visit(CRZ &crz) {}

void KernelToUnitaryVisitor::visit(CH &ch) {

  MatrixXcd fullMat = twoQubitGateExpand(CH_Mat, ch.bits()[0], ch.bits()[1]);

  m_circuitMat = fullMat * m_circuitMat;

}

void KernelToUnitaryVisitor::visit(S &s) {

  m_circuitMat = singleQubitGateExpand(S_Mat, s.bits()[0]) * m_circuitMat;

}

void KernelToUnitaryVisitor::visit(Sdg &sdg) {

  m_circuitMat = singleQubitGateExpand(Sdg_Mat, sdg.bits()[0]) * m_circuitMat;

}

void KernelToUnitaryVisitor::visit(T &t) {

  m_circuitMat = singleQubitGateExpand(T_Mat, t.bits()[0]) * m_circuitMat;

}

void KernelToUnitaryVisitor::visit(Tdg &tdg) {

  m_circuitMat = singleQubitGateExpand(Tdg_Mat, tdg.bits()[0]) * m_circuitMat;

}

void KernelToUnitaryVisitor::visit(CPhase &cphase) {xacc::error("Need to implement CPhase gate to matrix.");}



void KernelToUnitaryVisitor::visit(Measure &measure) {}

void KernelToUnitaryVisitor::visit(Identity &i) {}

void KernelToUnitaryVisitor::visit(U &u) {

  xacc::error("We don't support U3 gate to matrix.");

}

void KernelToUnitaryVisitor::visit(IfStmt &ifStmt) {}

// Identifiable Impl

MatrixXcd KernelToUnitaryVisitor::getMat() const { return m_circuitMat; }

MatrixXcd KernelToUnitaryVisitor::singleQubitGateExpand(MatrixXcd &in_gateMat,

                                                        size_t in_loc) const {

  if (m_nbQubit == 1) {

    return in_gateMat;

  }

  if (in_loc == 0) {

    auto dim = 1ULL << (m_nbQubit - 1);

    MatrixXcd rhs = MatrixXcd::Identity(dim, dim);

    return kroneckerProduct(in_gateMat, rhs);

  } else if (in_loc == m_nbQubit - 1) {

    auto dim = 1ULL << (m_nbQubit - 1);

    MatrixXcd lhs = MatrixXcd::Identity(dim, dim);

    return kroneckerProduct(lhs, in_gateMat);

  } else {

    auto leftDim = 1ULL << in_loc;

    auto rightDim = 1ULL << (m_nbQubit - 1 - in_loc);

    MatrixXcd lhs = MatrixXcd::Identity(leftDim, leftDim);

    MatrixXcd rhs = MatrixXcd::Identity(rightDim, rightDim);

    MatrixXcd lhsKron = kroneckerProduct(lhs, in_gateMat);

    return kroneckerProduct(lhsKron, rhs);

  }

}



MatrixXcd KernelToUnitaryVisitor::singleParametricQubitGateExpand(

    double in_var, Rot in_rotType, size_t in_bitLoc) const {

  static const std::complex<double> I_dual = I;

  MatrixXcd gateMat = [&]() {

    switch (in_rotType) {

      case Rot::X: {

        MatrixXcd rxMat =

            static_cast<std::complex<double>>(cos(in_var / 2)) *

                MatrixXcd::Identity(2, 2) -

            I_dual * static_cast<std::complex<double>>(sin(in_var / 2)) * X_Mat;

        return rxMat;

      }

      case Rot::Y: {

        MatrixXcd ryMat =

            static_cast<std::complex<double>>(cos(in_var / 2)) *

                MatrixXcd::Identity(2, 2) -

            I_dual * static_cast<std::complex<double>>(sin(in_var / 2)) * Y_Mat;

        return ryMat;

      }

      case Rot::Z: {

        MatrixXcd rzMat =

            static_cast<std::complex<double>>(cos(in_var / 2)) *

                MatrixXcd::Identity(2, 2) -

            I_dual * static_cast<std::complex<double>>(sin(in_var / 2)) * Z_Mat;

        return rzMat;

      }

      default:

        __builtin_unreachable();

    }

  }();

  return singleQubitGateExpand(gateMat, in_bitLoc);

}



MatrixXcd KernelToUnitaryVisitor::twoQubitGateExpand(MatrixXcd &in_gateMat,

                                                     size_t in_bit1,

                                                     size_t in_bit2) const {

  const auto matSize = 1ULL << m_nbQubit;

  MatrixXcd resultMat = MatrixXcd::Identity(matSize, matSize);

  // Qubit index list (0->(dim-1))

  std::vector<size_t> index_list(m_nbQubit);

  std::iota(index_list.begin(), index_list.end(), 0);

  const std::vector<size_t> idx{m_nbQubit - in_bit2 - 1,

                                m_nbQubit - in_bit1 - 1};

  for (size_t k = 0; k < matSize; ++k) {

    std::vector<std::complex<double>> oldColumn(matSize);

    for (size_t i = 0; i < matSize; ++i) {

      oldColumn[i] = resultMat(i, k);

    }



    for (size_t i = 0; i < matSize; ++i) {

      size_t local_i = 0;

      for (size_t l = 0; l < idx.size(); ++l) {

        local_i |= ((i >> idx[l]) & 1) << l;

      }



      std::complex<double> res = 0.0_i;

      for (size_t j = 0; j < (1ULL << idx.size()); ++j) {

        size_t locIdx = i;

        for (size_t l = 0; l < idx.size(); ++l) {

          if (((j >> l) & 1) != ((i >> idx[l]) & 1)) {

            locIdx ^= (1ULL << idx[l]);

          }

        }



        res += oldColumn[locIdx] * in_gateMat(local_i, j);

      }



      resultMat(i, k) = res;

    }

  }



  return resultMat;

}



}  // namespace qcor
 
 No newline at end of file