Merge pull request #104 from tnguyen-ornl/tnguyen/qcor-qs-callback

#pragma once

#include <thread>             

#include <mutex>              

#include <condition_variable> 

namespace qcor {

// Experimental optimization stepper to guide the Q# VQE loop.

class OptStepper {

public:

  OptStepper(const std::string &in_optName,

             xacc::HeterogeneousMap in_config = {}) {

    m_optimizer = qcor::createOptimizer(in_optName, std::move(in_config));

  }

  // Call by the main thread to update the  params and cost val...

  void update(const std::vector<double> &in_params, double in_costVal) {

    m_nextParamsAvail = false;

    if (m_dim == 0) {

      // First update:

      // Assuming the initial params are set by the caller,

      // i.e. this stepper doesn't control initial values.

      m_dim = in_params.size();

      m_optimizer->appendOption("initial-parameters", in_params);

      m_optFn = xacc::OptFunction(

          [&](const std::vector<double> &x, std::vector<double> &g) {

            return this->eval(x);

          },

          m_dim);

      m_currentCostVal = in_costVal;

      m_optParams = in_params;

      m_resultAvail = true;

      m_optThread = std::thread([&]() {

        auto result = m_optimizer->optimize(m_optFn);

        m_optDone = true;

      });

      return;

    }

    // A new iteration...

    m_optParams = in_params;

    m_currentCostVal = in_costVal;

    m_resultAvail = true;

    // Notify that the new result is avail...

    m_cv.notify_all();

  }

  // Fake evaluator function for the optimizer:

  // Run on the optimizer thread...

  double eval(const std::vector<double> &x) {

    if (!m_resultAvail) {

      m_optParams = x;

    }

    // Wait for the result to be available,

    // update() was called by the main thread.

    std::unique_lock<std::mutex> lck(m_mtx);

    m_cv.wait(lck, [this] { return m_resultAvail; });

    // This result has been processed...

    m_resultAvail = false;

    m_nextParamsAvail = true;

    return m_currentCostVal;

  }

  // Call by main thread...

  std::vector<double> getNextParams() {

    if (m_optDone) {

      // Return empty to signal that the optimizer has done.

      m_optThread.join();

      return {};

    }

    // Wait for the optimizer to give us a new set of parameters.

    // This should be fast, hence just poll.

    while (!m_nextParamsAvail) {

      std::this_thread::sleep_for(std::chrono::milliseconds(1));

    }

    return m_optParams;

  }

private:

  int m_dim = 0;

  std::shared_ptr<xacc::Optimizer> m_optimizer;

  std::vector<double> m_optParams;

  double m_currentCostVal;

  bool m_resultAvail = false;

  bool m_nextParamsAvail = false;

  bool m_optDone = false;

  // Optimizer thread

  std::mutex m_mtx;

  std::condition_variable m_cv;

  std::thread m_optThread;

  xacc::OptFunction m_optFn;

};

} // namespace qcor

namespace QCOR 

{

open QCOR.Intrinsic;

// Running VQE with an externally-provided optimization "stepper" interface:

// The stepper interface is "(current_params, energy) => (new_params)"

// Note: parameter initialization is done here.

// Returns the final energy.

operation DeuteronVqe(shots: Int, stepper : ((Double, Double[]) => Double[])) : Double {

    // Initial parameters

    let initial_params = [1.234];   

    mutable opt_params = initial_params;

    mutable energy_val = 0.0;

    use qubits = Qubit[2]

    {

        // Use repeat-until-success pattern:

        // when the optimization loop converges,

        // the stepper will return an empty param array.

        repeat {

            let xxExp = DeuteronXX(qubits, shots, opt_params[0]);

            let yyExp = DeuteronYY(qubits, shots, opt_params[0]);

            let z0_z1_exps = DeuteronZ0_Z1(qubits, shots, opt_params[0]);

            let z0Exp = z0_z1_exps[0];

            let z1Exp = z0_z1_exps[1];

            // Deuteron energy:

            // H = 5.907 - 2.1433 X0X1 - 2.1433 Y0Y1 + .21829 Z0 - 6.125 Z1

            set energy_val = 5.907 - 2.1433 * xxExp - 2.1433 * yyExp + 0.21829 * z0Exp - 6.125 * z1Exp;

            // Stepping...

            set opt_params = stepper(energy_val, opt_params);

        }

        until (Length(opt_params) == 0);

    }

    // Final energy:

    return energy_val;

}

// Base ansatz:

operation ansatz(qubits: Qubit[], theta: Double) : Unit {

    X(qubits[0]);

    Ry(theta, qubits[1]);

    CNOT(qubits[1], qubits[0]);

}

// This is for testing only:

// We should use a nicer implementation...

// XX term

operation DeuteronXX(qubits: Qubit[], shots: Int, theta: Double) : Double {

    mutable numParityOnes = 0;

    for shot in 1..shots {

        ansatz(qubits, theta);

        // Let's measure <X0X1>

        H(qubits[0]);

        H(qubits[1]);

        if M(qubits[0]) != M(qubits[1]) 

        {

            set numParityOnes += 1;

        }

        Reset(qubits[0]);

        Reset(qubits[1]);

    }

    let exp_val =  IntAsDouble(shots - numParityOnes)/IntAsDouble(shots) - IntAsDouble(numParityOnes)/IntAsDouble(shots);

    return exp_val;

}

// YY term

operation DeuteronYY(qubits: Qubit[], shots: Int, theta: Double) : Double {

    mutable numParityOnes = 0;

    for shot in 1..shots {

        ansatz(qubits, theta);

        // Let's measure <Y0Y1>

        Rx(1.57079632679, qubits[0]);

        Rx(1.57079632679, qubits[1]);

        if M(qubits[0]) != M(qubits[1]) 

        {

            set numParityOnes += 1;

        }

        Reset(qubits[0]);

        Reset(qubits[1]);

    }

    let exp_val =  IntAsDouble(shots - numParityOnes)/IntAsDouble(shots) - IntAsDouble(numParityOnes)/IntAsDouble(shots);

    return exp_val;

}

// Z0 and Z1 terms

operation DeuteronZ0_Z1(qubits: Qubit[], shots: Int, theta: Double) : Double[] {

    mutable numParityOnesZ0 = 0;

    mutable numParityOnesZ1 = 0;

    for shot in 1..shots {

        ansatz(qubits, theta);

        if M(qubits[0]) == One

        {

            set numParityOnesZ0 += 1;

        }

        if M(qubits[1]) == One

        {

            set numParityOnesZ1 += 1;

        }

        Reset(qubits[0]);

        Reset(qubits[1]);

    }

    let exp_val_z0 =  IntAsDouble(shots - numParityOnesZ0)/IntAsDouble(shots) - IntAsDouble(numParityOnesZ0)/IntAsDouble(shots);

    let exp_val_z1 =  IntAsDouble(shots - numParityOnesZ1)/IntAsDouble(shots) - IntAsDouble(numParityOnesZ1)/IntAsDouble(shots);

    return [exp_val_z0, exp_val_z1];

}

}

 No newline at end of file

#include <iostream> 

#include <vector>

#include "opt_stepper.hpp"

#include "qir-types-utils.hpp"

// Include the external QSharp function.

qcor_include_qsharp(QCOR__DeuteronVqe__body, double, int64_t shots, Callable* opt_stepper);

// Compile with:

// Include both the qsharp source and this driver file

// in the command line.

// $ qcor -qrt ftqc vqe_ansatz.qs vqe_driver.cpp

// Run with:

// $ ./a.out

int main() {

  // Create an optimizer:

  qcor::OptStepper qcorOptimizer("nlopt", {{"maxeval", 20}});

  using StepperFuncType =

      std::function<std::vector<double>(double, std::vector<double>)>;

  // Create an optimizer stepper as a lambda function to provide to Q#.

  StepperFuncType stepper_callable = [&](double in_costVal,

                                         std::vector<double> previous_params) {

    static size_t iterCount = 0;

    // Update the stepper with new data

    qcorOptimizer.update(previous_params, in_costVal);

    std::cout << "Iter " << iterCount++ << ": Energy(" << previous_params[0]

              << ") = " << in_costVal << "\n";

    // Returns a new set of params for Q# to try.

    return qcorOptimizer.getNextParams();

  };

  // Run the Q# Deuteron with the *nlopt* stepper provided as a callable.

  // Note: qsharp::createCallable will marshal the C++ function (lambda) to the

  // Q# Callable type.

  const int64_t nb_shots = 2048;

  const double final_energy = QCOR__DeuteronVqe__body(

      nb_shots, qir::createCallable(stepper_callable));

  std::cout << "Final energy = " << final_energy << "\n";

  return 0;

}

 No newline at end of file

namespace QCOR 

{

open QCOR.Intrinsic;

// Estimate energy value in a FTQC manner.

// Passing Array type b/w C++ driver and Q#

operation Ansatz(angles : Double[], shots: Int) : Double {

    mutable numParityOnes = 0;

    use qubits = Qubit[2]

    {

        for test in 1..shots {

            Rx(angles[0], qubits[0]);

            Rx(angles[1], qubits[1]);

            CNOT(qubits[1], qubits[0]);

            // Let's measure <X0X1>

            H(qubits[0]);

            H(qubits[1]);

            if M(qubits[0]) != M(qubits[1]) 

            {

                set numParityOnes += 1;

            }

            Reset(qubits[0]);

            Reset(qubits[1]);

        }

    }

    let res =  IntAsDouble(shots - numParityOnes)/IntAsDouble(shots) - IntAsDouble(numParityOnes)/IntAsDouble(shots);

    return res;

}

}

 No newline at end of file

#include <iostream>

#include <vector>

#include "qir-types-utils.hpp"

// Include the external QSharp function.

qcor_include_qsharp(QCOR__Ansatz__body, double, ::Array *, int64_t);

// Compile with:

// Include both the qsharp source and this driver file

// in the command line.

// $ qcor -qrt ftqc vqe_ansatz.qs vqe_driver.cpp

// Run with:

// $ ./a.out

int main() {

  const std::vector<double> angles{1.0, 2.0};

  const double exp_val_xx = QCOR__Ansatz__body(qir::toArray(angles), 1024);

  std::cout << "<X0X1> = " << exp_val_xx << "\n";

  return 0;

}

 No newline at end of file