aer_accelerator.cpp 33.6 KB
Newer Older
Mccaskey, Alex's avatar
Mccaskey, Alex committed
1
2
3
4
5
6
7
8
9
10
11
12
13
/*******************************************************************************
 * Copyright (c) 2019 UT-Battelle, LLC.
 * All rights reserved. This program and the accompanying materials
 * are made available under the terms of the Eclipse Public License v1.0
 * and Eclipse Distribution License v1.0 which accompanies this
 * distribution. The Eclipse Public License is available at
 * http://www.eclipse.org/legal/epl-v10.html and the Eclipse Distribution
 *License is available at https://eclipse.org/org/documents/edl-v10.php
 *
 * Contributors:
 *   Thien Nguyen - initial API and implementation
 *******************************************************************************/
#include "aer_accelerator.hpp"
14
#include "aer_noise_model.hpp"
Mccaskey, Alex's avatar
Mccaskey, Alex committed
15
16
17
18
19
20
21
22
23
24
25
26
27
28
#include "CommonGates.hpp"
#include "CountGatesOfTypeVisitor.hpp"
#include "InstructionIterator.hpp"

#include "controllers/controller.hpp"
#include "controllers/controller_execute.hpp"
#include "controllers/qasm_controller.hpp"
#include "controllers/statevector_controller.hpp"
#include "noise/readout_error.hpp"

#include "cppmicroservices/BundleActivator.h"
#include "cppmicroservices/BundleContext.h"
#include "cppmicroservices/ServiceProperties.h"
#include "xacc.hpp"
29
#include "xacc_service.hpp"
Mccaskey, Alex's avatar
Mccaskey, Alex committed
30
31

#include <bitset>
32
#include <optional>
33
#include "QObjGenerator.hpp"
34
#include "py-aer/aer_python_adapter.hpp"
Mccaskey, Alex's avatar
Mccaskey, Alex committed
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

namespace xacc {
namespace quantum {

#define IS_INTEGRAL(T)                                                         \
  typename std::enable_if<std::is_integral<T>::value>::type * = 0

template <class T>
std::string integral_to_binary_string(T byte, IS_INTEGRAL(T)) {
  std::bitset<sizeof(T) * 8> bs(byte);
  return bs.to_string();
}
std::string hex_string_to_binary_string(std::string hex) {
  return integral_to_binary_string((int)strtol(hex.c_str(), NULL, 0));
}

51
HeterogeneousMap AerAccelerator::getProperties() {
52
53
  auto props = physical_backend_properties;
  // Insert 'shots' data
54
55
56
  if (m_simtype == "qasm") {
    props.insert("shots", m_shots);
  }
57
  return props;
58
59
}

Mccaskey, Alex's avatar
Mccaskey, Alex committed
60
61
void AerAccelerator::initialize(const HeterogeneousMap &params) {

62
  m_options = params;
Mccaskey, Alex's avatar
Mccaskey, Alex committed
63
64
  noise_model.clear();
  m_simtype = "qasm";
65
  connectivity.clear();
Mccaskey, Alex's avatar
Mccaskey, Alex committed
66
67
68
69
70
71
72

  xacc_to_qobj = xacc::getCompiler("qobj");
  if (params.keyExists<int>("shots")) {
    m_shots = params.get<int>("shots");
  }
  if (params.stringExists("sim-type")) {
    if (!xacc::container::contains(
73
            std::vector<std::string>{"qasm", "statevector", "pulse", "density_matrix"},
Mccaskey, Alex's avatar
Mccaskey, Alex committed
74
75
76
77
78
79
80
81
82
83
            params.getString("sim-type"))) {
      xacc::warning("[Aer] warning, invalid sim-type (" +
                    params.getString("sim-type") +
                    "), must be qasm or statevector.");
    } else {
      m_simtype = params.getString("sim-type");
    }
  }

  if (params.stringExists("backend")) {
84
85
86
87
88
89
    auto ibm_noise_model = xacc::getService<NoiseModel>("IBM");
    ibm_noise_model->initialize(params);
    auto json_str = ibm_noise_model->toJson();
    noise_model = nlohmann::json::parse(json_str);
    auto ibm = xacc::getAccelerator("ibm:" + params.getString("backend"));
    physical_backend_properties = ibm->getProperties();
90
91
92
93
    if (m_simtype == "pulse") {
      // If pulse mode, must contribute pulse cmd-def.
      ibm->contributeInstructions();
    }
94
95
    // Set connectivity based on the backend:
    connectivity = ibm->getConnectivity();
96
  } else if (params.stringExists("noise-model")) {
97
98
99
100
101
102
103
104
105
106
    std::string noise_model_str = params.getString("noise-model");
    // Check if this is a file name
    std::ifstream test(noise_model_str);
    if (test) {
      std::string str((std::istreambuf_iterator<char>(test)),
                      std::istreambuf_iterator<char>());
      noise_model = nlohmann::json::parse(str);
    } else {
      noise_model = nlohmann::json::parse(params.getString("noise-model"));
    }
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
    
    // need this for ro error decorator
    std::vector<double> p01s, p10s;
    auto errors = noise_model["errors"];
    for (auto error : errors) {
      if (error["operations"].get<std::vector<std::string>>() ==
          std::vector<std::string>{"measure"}) {
        auto probs =
            error["probabilities"].get<std::vector<std::vector<double>>>();
        auto meas1prep0 = probs[0][1];
        auto meas0prep1 = probs[1][0];
        p01s.push_back(meas0prep1);
        p10s.push_back(meas1prep0);
      }
    }
    physical_backend_properties.insert("p01s", p01s);
    physical_backend_properties.insert("p10s", p10s);

Mccaskey, Alex's avatar
Mccaskey, Alex committed
125
  }
126
  initialized = true;
Mccaskey, Alex's avatar
Mccaskey, Alex committed
127
}
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
double AerAccelerator::calcExpectationValueZ(
    const std::vector<std::pair<double, double>> &in_stateVec,
    const std::vector<std::size_t> &in_bits) {
  const auto hasEvenParity =
      [](size_t x, const std::vector<size_t> &in_qubitIndices) -> bool {
    size_t count = 0;
    for (const auto &bitIdx : in_qubitIndices) {
      if (x & (1ULL << bitIdx)) {
        count++;
      }
    }
    return (count % 2) == 0;
  };

  double result = 0.0;
  for (uint64_t i = 0; i < in_stateVec.size(); ++i) {
    result += (hasEvenParity(i, in_bits) ? 1.0 : -1.0) *
              std::norm(std::complex<double>(in_stateVec[i].first,
                                             in_stateVec[i].second));
  }

  return result;
}
Mccaskey, Alex's avatar
Mccaskey, Alex committed
151

152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
double AerAccelerator::calcExpectationValueZFromDensityMatrix(
    const std::vector<std::vector<std::pair<double, double>>> &in_densityMat,
    const std::vector<std::size_t> &in_bits) {
  const auto hasEvenParity =
      [](size_t x, const std::vector<size_t> &in_qubitIndices) -> bool {
    size_t count = 0;
    for (const auto &bitIdx : in_qubitIndices) {
      if (x & (1ULL << bitIdx)) {
        count++;
      }
    }
    return (count % 2) == 0;
  };

  double result = 0.0;
  for (uint64_t i = 0; i < in_densityMat.size(); ++i) {
    const auto &diag_elem = in_densityMat[i][i];
    // The diag. elements of the DM should be real.
    assert(std::abs(diag_elem.second) < 1e-3);
    // When using DM elements, don't need the square-norm.
    result += ((hasEvenParity(i, in_bits) ? 1.0 : -1.0) * diag_elem.first);
  }

  return result;
}

Mccaskey, Alex's avatar
Mccaskey, Alex committed
178
179
180
181
182
183
184
185
186
187
void AerAccelerator::execute(
    std::shared_ptr<AcceleratorBuffer> buffer,
    const std::shared_ptr<CompositeInstruction> program) {

  if (m_simtype == "qasm") {
    auto qobj_str = xacc_to_qobj->translate(program);

    nlohmann::json j = nlohmann::json::parse(qobj_str)["qObject"];
    j["config"]["shots"] = m_shots;
    j["config"]["noise_model"] = noise_model;
188
189
190

    // xacc::set_verbose(true);
    // xacc::info("Shots Qobj:\n" + j.dump(2));
Mccaskey, Alex's avatar
Mccaskey, Alex committed
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
    auto results_json = nlohmann::json::parse(
        AER::controller_execute_json<AER::Simulator::QasmController>(j.dump()));

    auto results = *results_json["results"].begin();

    auto counts = results["data"]["counts"].get<std::map<std::string, int>>();
    for (auto &kv : counts) {
      std::string hexStr = kv.first;
      int nOccurrences = kv.second;
      auto bitStr = hex_string_to_binary_string(hexStr);

      // Process bitStr to be an n-Measure string in msb
      std::string actual = "";
      CountGatesOfTypeVisitor<Measure> cc(program);
      int nMeasures = cc.countGates();
      for (int i = 0; i < nMeasures; i++)
        actual += "0";
      for (int i = 0; i < nMeasures; i++)
        actual[actual.length() - 1 - i] = bitStr[bitStr.length() - i - 1];

      buffer->appendMeasurement(actual, nOccurrences);
    }
213
214
215
216
217
218
219
220
221
222
  } else if (m_simtype == "pulse") {
    // Get the correct QObject Generator
    auto qobjGen = xacc::getService<QObjGenerator>("pulse");
    auto chosenBackend = nlohmann::json::parse(physical_backend_properties.getString("config-json"));
    // Unused
    const std::string getBackendPropsResponse = "{}";
    auto defaults_response = nlohmann::json::parse(physical_backend_properties.getString("defaults-json"));
    auto ibmPulseAssembler = xacc::getService<IRTransformation>("ibm-pulse");
    auto kernel = xacc::ir::asComposite(program->clone());
    
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
    // Remove measures, add measure all (at the end)
    kernel->clear();
    InstructionIterator iter(program);
    std::vector<size_t> measured_bits;
    while (iter.hasNext()) {
      auto next = iter.next();
      if (!next->isComposite() && next->name() != "Measure") {
        kernel->addInstruction(next);
      } else if (next->name() == "Measure") {
        measured_bits.push_back(next->bits()[0]);
      }
    }

    auto provider = xacc::getIRProvider("quantum");
    for (size_t qId = 0; qId < buffer->size(); ++qId) {
      kernel->addInstruction(provider->createInstruction("Measure", qId));
    }

    // Assemble pulse composite from the input kernel.
    ibmPulseAssembler->apply(kernel, nullptr);    
243
    // Generate the QObject JSON
244
    auto qobjJsonStr = qobjGen->getQObjJsonStr({kernel}, m_shots, chosenBackend,
245
246
                                           getBackendPropsResponse,
                                           connectivity, defaults_response);
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
    xacc::info("Qobj:\n" + qobjJsonStr);
    auto hamiltonianJson = chosenBackend["hamiltonian"];
    // Remove unrelated fields (could contain problematic characters)
    hamiltonianJson.erase("description");
    hamiltonianJson.erase("h_latex"); 
    xacc::info("Hamiltonian Json:\n" + hamiltonianJson.dump());
    const auto dt = chosenBackend["dt"].get<double>();
    const auto qubitFreqEst = defaults_response["qubit_freq_est"].get<std::vector<double>>();
    const auto uLoFreqs = chosenBackend["u_channel_lo"];
    std::vector<int> uLoRefs;
    for (auto loIter = uLoFreqs.begin(); loIter != uLoFreqs.end(); ++loIter) {
      auto uLoConfig = *((*loIter).begin());
      uLoRefs.emplace_back(uLoConfig["q"].get<int>());
    }
    // Run the simulation via Python
    const std::string resultJson =
        xacc::aer::runPulseSim(hamiltonianJson.dump(), dt, qubitFreqEst, uLoRefs, qobjJsonStr);
264
265
    auto result_json = nlohmann::json::parse(resultJson);
    auto count_json = result_json["counts"].get<std::map<std::string, int>>();
266
267
268
269
270
271
272
273
274
275
276
277
    for (const auto &[hexStr, nOccurrences] : count_json) {
      auto bitStr = hex_string_to_binary_string(hexStr);
      // Process bitStr to be an n-Measure string in msb
      std::string actual = "";
      const auto nMeasures = measured_bits.size();
      for (int i = 0; i < nMeasures; i++) {
        actual += "0";
      }

      for (int i = 0; i < nMeasures; i++) {
        actual[actual.length() - 1 - i] = bitStr[bitStr.length() - measured_bits[i] - 1];
      }
278

279
280
      buffer->appendMeasurement(actual, nOccurrences);
    }
281
282
283

    auto state_vector = result_json["statevector"].get<std::vector<std::pair<double, double>>>();
    buffer->addExtraInfo("state", state_vector);
284
285
286
287
288
289
290
291
292
293
294
295
296
297
  } else if (m_simtype == "density_matrix") {
    // remove all measures, don't need them
    auto tmp = xacc::ir::asComposite(program->clone());
    tmp->clear();
    std::vector<std::size_t> measured_bits;
    InstructionIterator iter(program);
    while (iter.hasNext()) {
      auto next = iter.next();
      if (!next->isComposite() && next->name() != "Measure") {
        tmp->addInstruction(next);
      } else if (next->name() == "Measure") {
        measured_bits.push_back(next->bits()[0]);
      }
    }
298
299
300
    // In "density_matrix" simulation mode, always include Id gates
    // if they are explicitly added to the Composite.
    auto qobj_str = xacc_to_qobj->translate(tmp, {{"skip-id-gates", false}});
301
    nlohmann::json j = nlohmann::json::parse(qobj_str)["qObject"];
302
    xacc::info("Qobj:\n" + j.dump());
303
304
305
306
307
308
309
310
    j["config"]["noise_model"] = noise_model;
    j["config"]["method"] = "density_matrix";
    auto snapshotInst = nlohmann::json::object();
    snapshotInst["label"] = "dm_snapshot";
    snapshotInst["name"] = "snapshot";
    snapshotInst["snapshot_type"] = "density_matrix";
    auto& exprJson = *(j["experiments"].begin());
    exprJson["instructions"].push_back(snapshotInst);
311
    // std::cout << "Qobj:\n" << j.dump(2);
312
313
    auto results_json = nlohmann::json::parse(
        AER::controller_execute_json<AER::Simulator::QasmController>(j.dump()));
314
    // std::cout << "Result:\n" << results_json.dump() << "\n";
315
316
317
318
319
320
321
322
323
324
    auto results = *results_json["results"].begin();
    auto dm_mat =
        (*(results["data"]["snapshots"]["density_matrix"]["dm_snapshot"]
               .begin()))["value"]
            .get<std::vector<std::vector<std::pair<double, double>>>>();
    std::vector<std::pair<double, double>> flattenDm;
    for (const auto &row : dm_mat) {
      flattenDm.insert(flattenDm.end(), row.begin(), row.end());
    }
    buffer->addExtraInfo("density_matrix", flattenDm);
325
326
    auto exp_val = calcExpectationValueZFromDensityMatrix(dm_mat, measured_bits);
    buffer->addExtraInfo("exp-val-z", exp_val);
327
328
  }
  else {
Mccaskey, Alex's avatar
Mccaskey, Alex committed
329
330
331
332
    // statevector
    // remove all measures, don't need them
    auto tmp = xacc::ir::asComposite(program->clone());
    tmp->clear();
333
    std::vector<std::size_t> measured_bits;
Mccaskey, Alex's avatar
Mccaskey, Alex committed
334
335
336
337
338
    InstructionIterator iter(program);
    while (iter.hasNext()) {
      auto next = iter.next();
      if (!next->isComposite() && next->name() != "Measure") {
        tmp->addInstruction(next);
339
340
      } else if (next->name() == "Measure") {
        measured_bits.push_back(next->bits()[0]);
Mccaskey, Alex's avatar
Mccaskey, Alex committed
341
342
343
344
345
      }
    }
    auto qobj_str = xacc_to_qobj->translate(tmp);

    nlohmann::json j = nlohmann::json::parse(qobj_str)["qObject"];
346
347
    j["config"]["noise_model"] = noise_model;
    // xacc::info("StateVec Qobj:\n" + j.dump(2));
Mccaskey, Alex's avatar
Mccaskey, Alex committed
348
349
350
351
352

    auto results_json = nlohmann::json::parse(
        AER::controller_execute_json<AER::Simulator::StatevectorController>(
            j.dump()));

353
354
355
356
    if (results_json["status"].get<std::string>().find("ERROR") !=
        std::string::npos) {
      std::cout << results_json["status"].get<std::string>() << "\n";
      xacc::error("Aer Error: " + results_json["status"].get<std::string>());
357
358
    }

Mccaskey, Alex's avatar
Mccaskey, Alex committed
359
360
361
362
363
364
365
    auto results = *results_json["results"].begin();

    auto state_vector = results["data"]["statevector"]
                            .get<std::vector<std::pair<double, double>>>();

    buffer->addExtraInfo("state", state_vector);

366
367
    auto exp_val = calcExpectationValueZ(state_vector, measured_bits);
    buffer->addExtraInfo("exp-val-z", exp_val);
Mccaskey, Alex's avatar
Mccaskey, Alex committed
368
369
370
371
372
373
374
375
376
377
378
379
380
381
  }
}

void AerAccelerator::execute(
    std::shared_ptr<AcceleratorBuffer> buffer,
    const std::vector<std::shared_ptr<CompositeInstruction>>
        compositeInstructions) {
  for (auto &f : compositeInstructions) {
    auto tmpBuffer =
        std::make_shared<xacc::AcceleratorBuffer>(f->name(), buffer->size());
    execute(tmpBuffer, f);
    buffer->appendChild(f->name(), tmpBuffer);
  }
}
382

383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
void AerAccelerator::apply(std::shared_ptr<AcceleratorBuffer> buffer,
                           std::shared_ptr<Instruction> inst) {
  static AER::Statevector::State<QV::QubitVector<double>> stateVec;
  static auto provider = xacc::getIRProvider("quantum");
  static AER::RngEngine rng(time(NULL));

  if (!noiseModelObj) {
    noiseModelObj = std::make_shared<AER::Noise::NoiseModel>(noise_model);
    stateVec.initialize_qreg(buffer->size());
  }
  if (inst->isComposite() || inst->isAnalog()) {
    xacc::error("Only gates are allowed.");
  }

  auto tempComp = provider->createComposite("tmp");
  tempComp->addInstruction(inst);
  auto qobj_str = xacc_to_qobj->translate(tempComp);
  auto qObjJson = nlohmann::json::parse(qobj_str)["qObject"];
  AER::Qobj qobj(qObjJson);
  assert(qobj.circuits.size() == 1);
  auto circ = qobj.circuits[0];
404

405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
  // Output data container
  AER::ExperimentData data;
  data.add_metadata("method", stateVec.name());
  data.add_metadata("measure_sampling", false);
  auto noiseCirc = noiseModelObj->sample_noise(circ, rng);
  stateVec.initialize_creg(circ.num_memory, circ.num_registers);
  stateVec.apply_ops(noiseCirc.ops, data, rng);
  stateVec.add_creg_to_data(data);
  // If it was a Measure op:
  if (inst->name() == "Measure") {
    xacc::info("Experiment data: \n" + data.json().dump() + "\n");
    auto countData = data.json()["counts"].get<std::map<std::string, int>>();
    // In this mode, we only measure 1 qubit for 1 shot
    assert(countData.size() == 1);
    int zeroCount = countData.find("0x0") == countData.end() ? 0 : 1;
    int oneCount = countData.find("0x1") == countData.end() ? 0 : 1;
    assert((zeroCount + oneCount) == 1);
    buffer->measure(inst->bits()[0], (oneCount > 0 ? 1 : 0));
  }
}

426
427
428
429
430
431
432
void IbmqNoiseModel::initialize(const HeterogeneousMap &params) {
  m_nbQubits = 0;
  m_qubitT1.clear();
  m_qubitT2.clear();
  m_gateErrors.clear();
  m_gateDurations.clear();
  m_roErrors.clear();
433
  m_connectivity.clear();
434
435
436
437
438
439
440
  // Note: we support both remote backend JSON and cache JSON string.
  // So that we can test this with offline JSON.
  if (params.stringExists("backend") || params.stringExists("backend-json")) {
    const auto backendJsonStr = [&]() {
      if (params.stringExists("backend")) {
        auto ibm = xacc::getAccelerator("ibm:" + params.getString("backend"));
        auto props = ibm->getProperties().get<std::string>("total-json");
441
        m_connectivity = ibm->getConnectivity();
442
443
444
445
446
447
        return props;
      }
      return params.getString("backend-json");
    }();

    auto backEndJson = nlohmann::json::parse(backendJsonStr);
448
449
    // Cache the backend properties JSON.
    m_backendPropertiesJson = backendJsonStr;
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
    // Parse qubit data:
    auto qubitsData = backEndJson["qubits"];
    size_t nbQubit = 0;
    for (auto qubitIter = qubitsData.begin(); qubitIter != qubitsData.end();
         ++qubitIter) {
      std::optional<double> meas0Prep1, meas1Prep0;
      // Each qubit contains an array of properties.
      for (auto probIter = qubitIter->begin(); probIter != qubitIter->end();
           ++probIter) {
        const auto probObj = *probIter;
        const std::string probName = probObj["name"].get<std::string>();
        const double probVal = probObj["value"].get<double>();
        const std::string unit = probObj["unit"].get<std::string>();
        if (probName == "T1") {
          assert(unit == "µs" || unit == "us" || unit == "ns");
465
466
          m_qubitT1.emplace_back(
              (unit == "µs" || unit == "us") ? 1000.0 * probVal : probVal);
467
468
469
470
        }

        if (probName == "T2") {
          assert(unit == "µs" || unit == "us" || unit == "ns");
471
472
          m_qubitT2.emplace_back(
              (unit == "µs" || unit == "us") ? 1000.0 * probVal : probVal);
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
        }

        if (probName == "prob_meas0_prep1") {
          assert(unit.empty());
          meas0Prep1 = probVal;
        }

        if (probName == "prob_meas1_prep0") {
          assert(unit.empty());
          meas1Prep0 = probVal;
        }
      }
      assert(meas0Prep1.has_value() && meas1Prep0.has_value());
      m_roErrors.emplace_back(std::make_pair(*meas0Prep1, *meas1Prep0));

      nbQubit++;
    }
    m_nbQubits = nbQubit;
    assert(m_qubitT1.size() == m_nbQubits);
    assert(m_qubitT2.size() == m_nbQubits);
    assert(m_roErrors.size() == m_nbQubits);

    // Parse gate data:
    auto gateData = backEndJson["gates"];
    for (auto gateIter = gateData.begin(); gateIter != gateData.end();
         ++gateIter) {
      auto gateObj = *gateIter;
      const std::string gateName = gateObj["name"].get<std::string>();
      auto gateParams = gateObj["parameters"];
      for (auto it = gateParams.begin(); it != gateParams.end(); ++it) {
        auto paramObj = *it;
        const std::string paramName = paramObj["name"].get<std::string>();
        if (paramName == "gate_length") {
          const std::string unit = paramObj["unit"].get<std::string>();
          assert(unit == "µs" || unit == "us" || unit == "ns");
          const double gateLength =
509
510
511
              (unit == "µs" || unit == "us")
                  ? 1000.0 * paramObj["value"].get<double>()
                  : paramObj["value"].get<double>();
512
513
          const bool insertOk =
              m_gateDurations.insert({gateName, gateLength}).second;
514
515
516
517
518
519
520
521
522
523
524
525
526
527
          if (gateName.rfind("sx", 0) == 0) {
            // This device properties were specified in the { rz, sx, cx } basis.
            // We add equivalent data for { u1, u2, u3 } basis set as well.
            // For every sx gate, add equiv. data for u1, u2, u3 gates.
            // Note: some of our noisy simulators don't split a gate into multiple sx gates,
            // hence, they cannot simulate noise at that level.
            const std::string qubitOperandSuffix = gateName.substr(2);
            const std::string u1GateName = "u1" + qubitOperandSuffix;
            const std::string u2GateName = "u2" + qubitOperandSuffix;
            const std::string u3GateName = "u3" + qubitOperandSuffix;
            m_gateDurations.insert({u1GateName, 0.0});
            m_gateDurations.insert({u2GateName, gateLength});
            m_gateDurations.insert({u3GateName, 2.0 * gateLength});
          }
528
529
530
531
532
533
534
535
536
537
538
          // Must not contain duplicates.
          assert(insertOk);
        }

        if (paramName == "gate_error") {
          assert(paramObj["unit"].get<std::string>().empty());
          const bool insertOk =
              m_gateErrors.insert({gateName, paramObj["value"].get<double>()})
                  .second;
          // Must not contain duplicates.
          assert(insertOk);
539
540
541
542
543
544
545
546
547
          if (gateName.rfind("sx", 0) == 0) {
            const std::string qubitOperandSuffix = gateName.substr(2);
            const std::string u1GateName = "u1" + qubitOperandSuffix;
            const std::string u2GateName = "u2" + qubitOperandSuffix;
            const std::string u3GateName = "u3" + qubitOperandSuffix;
            m_gateErrors.insert({u1GateName, 0.0});
            m_gateErrors.insert({u2GateName, paramObj["value"].get<double>()});
            m_gateErrors.insert({u3GateName, 2.0 * paramObj["value"].get<double>()});
          }
548
549
550
551
552
553
        }
      }
    }
  }
}

554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
std::string
IbmqNoiseModel::getUniversalGateEquiv(xacc::quantum::Gate &in_gate) const {
  if (in_gate.bits().size() == 1 && in_gate.name() != "Measure") {
    // Note: rotation around Z is a noiseless *u1* operation;
    // *u2* operations are those that requires a half-length rotation;
    // *u3* operations are those that requires a full-length rotation.
    static const std::unordered_map<std::string, std::string>
        SINGLE_QUBIT_GATE_MAP{{"X", "u3"},   {"Y", "u3"},  {"Z", "u1"},
                              {"H", "u2"},   {"U", "u3"},  {"T", "u1"},
                              {"Tdg", "u1"}, {"S", "u1"},  {"Sdg", "u1"},
                              {"Rz", "u1"},  {"Rx", "u3"}, {"Ry", "u3"}};
    const auto iter = SINGLE_QUBIT_GATE_MAP.find(in_gate.name());
    // If cannot find the gate, just treat that as a noiseless u1 op.
    const std::string universalGateName =
        (iter == SINGLE_QUBIT_GATE_MAP.end()) ? "u1" : iter->second;
569
    return universalGateName + std::to_string(in_gate.bits()[0]);
570
571
572
573
574
575
576
  }

  if (in_gate.bits().size() == 2) {
    return "cx" + std::to_string(in_gate.bits()[0]) + "_" +
           std::to_string(in_gate.bits()[1]);
  }

577
  return "id" + std::to_string(in_gate.bits()[0]);
578
579
}

580
581
582
583
// Return gate time, T1, T2
std::tuple<double, double, double>
IbmqNoiseModel::relaxationParams(xacc::quantum::Gate &in_gate,
                                 size_t in_qubitIdx) const {
584
585
586
587
  const std::string universalGateName = getUniversalGateEquiv(in_gate);
  const auto gateDurationIter = m_gateDurations.find(universalGateName);
  assert(gateDurationIter != m_gateDurations.end());
  const double gateDuration = gateDurationIter->second;
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
  const double qubitT1 = m_qubitT1[in_qubitIdx];
  const double qubitT2 = m_qubitT2[in_qubitIdx];
  return std::make_tuple(gateDuration, qubitT1, qubitT2);
}

std::vector<std::vector<std::complex<double>>>
IbmqNoiseModel::thermalRelaxationChoiMat(double in_gateTime, double in_T1,
                                         double in_T2) const {
  const double rate1 = 1.0 / in_T1;
  const double pReset = 1.0 - std::exp(-in_gateTime * rate1);

  const double rate2 = 1.0 / in_T2;
  const double expT2 = std::exp(-in_gateTime * rate2);
  const double p0 = 1.0;
  const double p1 = 0.0;
  return {{1 - p1 * pReset, 0, 0, expT2},
          {0, p1 * pReset, 0, 0},
          {0, 0, p0 * pReset, 0},
          {expT2, 0, 0, 1 - p0 * pReset}};
607
608
609
610
}

std::vector<double> IbmqNoiseModel::calculateDepolarizing(
    xacc::quantum::Gate &in_gate,
611
612
    const std::vector<std::vector<std::complex<double>>> &in_relaxationError)
    const {
613
614
615
616
617
618
  //  Compute the depolarizing channel error parameter in the
  //  presence of T1/T2 thermal relaxation.
  //  Hence we have that the depolarizing error probability
  //  for the composed depolarization channel is
  //  p = dim * (F(E_relax) - F) / (dim * F(E_relax) - 1)
  const double averageThermalError =
619
      1.0 - averageGateFidelity(in_gate, in_relaxationError);
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
  const std::string universalGateName = getUniversalGateEquiv(in_gate);
  // Retrieve the error rate:
  const auto gateErrorIter = m_gateErrors.find(universalGateName);
  const double gateError =
      (gateErrorIter == m_gateErrors.end()) ? 0.0 : gateErrorIter->second;
  // If the backend gate error (estimated by randomized benchmarking) is more
  // than thermal relaxation error. We need to add depolarization to simulate
  // the total gate error.
  if (gateError > averageThermalError) {
    // Model gate error entirely as depolarizing error
    const double depolError = 2 * (gateError - averageThermalError) /
                              (2 * (1 - averageThermalError) - 1);
    return {depolError};
  }
  return {0.0};
}

double IbmqNoiseModel::averageGateFidelity(
    xacc::quantum::Gate &in_gate,
639
640
    const std::vector<std::vector<std::complex<double>>> &in_relaxationError)
    const {
641
  // We only handle single-qubit gates for now.
642
  if (in_gate.bits().size() == 1 && in_relaxationError.size() == 4) {
643
644
    // Note: this is based on the simplified amplitude-damping channel
    const double processFidelity =
645
646
        (in_relaxationError[0][0].real() + in_relaxationError[0][3].real() +
         in_relaxationError[3][0].real() + in_relaxationError[3][3].real()) /
647
        4.0;
648

649
650
651
652
653
654
655
    assert(processFidelity <= 1.0);
    const double averageFidelity = (4.0 * processFidelity + 1.0) / 5.0;
    return averageFidelity;
  }
  return 1.0;
}

656
657
658
659
660
661
double IbmqNoiseModel::gateErrorProb(xacc::quantum::Gate &gate) const {
  const auto universalGateName = getUniversalGateEquiv(gate);
  const auto gateErrorIter = m_gateErrors.find(universalGateName);
  return (gateErrorIter == m_gateErrors.end()) ? 0.0 : gateErrorIter->second;
}

662
663
664
std::vector<NoiseChannelKraus>
IbmqNoiseModel::getNoiseChannels(xacc::quantum::Gate &gate) const {
  std::vector<NoiseChannelKraus> krausOps;
665
  const auto noiseUtils = xacc::getService<NoiseModelUtils>("default");
666
  if (gate.bits().size() == 1 && gate.name() != "Measure") {
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
667
    // Amplitude damping + dephasing
668
669
670
671
    const auto [gateDuration, qubitT1, qubitT2] =
        relaxationParams(gate, gate.bits()[0]);
    const auto relaxationError =
        thermalRelaxationChoiMat(gateDuration, qubitT1, qubitT2);
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
672
673

    // Depolarization
674
    const auto dpAmpl = calculateDepolarizing(gate, relaxationError);
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
675
676
    if (!dpAmpl.empty()) {
      const double probDP = dpAmpl[0];
677
678
679
680
681
      const std::vector<std::vector<std::complex<double>>> depolErrorChoi{
          {1.0 - probDP / 2.0, 0., 0., 1.0 - probDP},
          {0., probDP / 2.0, 0., 0.},
          {0., 0., probDP / 2.0, 0.},
          {1.0 - probDP, 0., 0., 1.0 - probDP / 2.0}};
682
      const auto noiseUtils = xacc::getService<NoiseModelUtils>("default");
683
      krausOps.emplace_back(NoiseChannelKraus(gate.bits(), noiseUtils->choiToKraus(depolErrorChoi), KrausMatBitOrder::MSB));
684
685
    }
  }
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
686
687
  // For two-qubit gates, we currently only support
  // amplitude damping on both qubits (scaled by gate time).
688
  // We don't have ability to handle multi-qubit depolarization
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
689
690
  // in TNQVM yet.
  if (gate.bits().size() == 2) {
691
692
693
694
695
696

    for (const auto &qubitIdx : gate.bits()) {
      const auto [gateDuration, qubitT1, qubitT2] =
          relaxationParams(gate, qubitIdx);
      const auto relaxationError =
          thermalRelaxationChoiMat(gateDuration, qubitT1, qubitT2);
697
      krausOps.emplace_back(NoiseChannelKraus({qubitIdx}, noiseUtils->choiToKraus(relaxationError), KrausMatBitOrder::MSB));
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
698
699
    }
  }
700
701
702
703
  return krausOps;
}

std::string IbmqNoiseModel::toJson() const {
704
705
706
707
708
709
710
711
712
713
  // First, check if Qiskit is available,
  // try to use Qiskit util if possible.
  const std::string noiseModelJson =
      xacc::aer::noiseModelFromBackendProperties(m_backendPropertiesJson);
  if (!noiseModelJson.empty()) {
    xacc::info("Qiskit generated noise model:\n" + noiseModelJson);
    return noiseModelJson;
  }
  
  // No Qiskit, create the noise model by ourselves from the backend properties.
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
  // Aer noise model Json
  nlohmann::json noiseModel;
  std::vector<nlohmann::json> noiseElements;
  // Adds RO errors:
  const auto roErrors = readoutErrors();
  for (size_t qIdx = 0; qIdx < roErrors.size(); ++qIdx) {
    const auto &[meas0Prep1, meas1Prep0] = roErrors[qIdx];
    const std::vector<std::vector<double>> probs{{1 - meas1Prep0, meas1Prep0},
                                                 {meas0Prep1, 1 - meas0Prep1}};
    nlohmann::json element;
    element["type"] = "roerror";
    element["operations"] = std::vector<std::string>{"measure"};
    element["probabilities"] = probs;
    element["gate_qubits"] = std::vector<std::vector<std::size_t>>{{qIdx}};
    noiseElements.push_back(element);
  }

731
  const auto noiseUtils = xacc::getService<NoiseModelUtils>("default");
732
  // Add Kraus noise:
733
  // (1) Single-qubit gate noise:
734
735
  // Note: we must add noise ops for u2, u3, and cx gates:
  for (size_t qIdx = 0; qIdx < roErrors.size(); ++qIdx) {
736
737
738
    // For mapping purposes:
    // U2 == Hadamard gate
    // U3 == X gate
739
740
    Hadamard gateU2(qIdx);
    X gateU3(qIdx);
741
742
    const std::unordered_map<std::string, xacc::quantum::Gate *> gateMap{
        {"u2", &gateU2}, {"u3", &gateU3}};
743
744

    for (const auto &[gateName, gate] : gateMap) {
745
      const auto errorChannels = getNoiseChannels(*gate);
746
747
748
749
      nlohmann::json element;
      element["type"] = "qerror";
      element["operations"] = std::vector<std::string>{gateName};
      element["gate_qubits"] = std::vector<std::vector<std::size_t>>{{qIdx}};
750
      std::vector<nlohmann::json> krausOps;
751
      for (const auto &error : errorChannels) {
752
        const auto krausOpMats = error.mats;
753
754
755
756
757
758
        nlohmann::json instruction;
        instruction["name"] = "kraus";
        instruction["qubits"] = std::vector<std::size_t>{0};
        instruction["params"] = krausOpMats;
        krausOps.emplace_back(instruction);
      }
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
      element["instructions"] =
          std::vector<std::vector<nlohmann::json>>{krausOps};
      element["probabilities"] = std::vector<double>{1.0};
      noiseElements.push_back(element);
    }
  }

  // (2) Two-qubit gate noise:
  for (const auto &[qubit1, qubit2] : m_connectivity) {
    // We need to add noise for both CNOT directions
    // Note: the duration of them can be different,
    // hence the noise channels.
    CNOT cx1(std::vector<std::size_t>{(size_t)qubit1, (size_t)qubit2});
    CNOT cx2(std::vector<std::size_t>{(size_t)qubit2, (size_t)qubit1});
    const std::vector<xacc::quantum::Gate *> cxGates{&cx1, &cx2};
    for (const auto &cx : cxGates) {
775
      const auto errorChannels = getNoiseChannels(*cx);
776
777
778
779
780
781
782
783
784
785
      assert(errorChannels.size() == 2);
      nlohmann::json element;
      element["type"] = "qerror";
      element["operations"] =
          std::vector<std::string>{getUniversalGateEquiv(*cx)};
      element["gate_qubits"] =
          std::vector<std::vector<std::size_t>>{cx->bits()};
      std::vector<nlohmann::json> krausOps;
      size_t noiseBitIdx = 0;
      for (const auto &error : errorChannels) {
786
        const auto krausOpMats = error.mats;
787
788
789
790
791
792
793
794
795
796
        nlohmann::json instruction;
        instruction["name"] = "kraus";
        instruction["qubits"] = std::vector<std::size_t>{noiseBitIdx++};
        instruction["params"] = krausOpMats;
        krausOps.emplace_back(instruction);
      }
      assert(krausOps.size() == 2);
      element["instructions"] =
          std::vector<std::vector<nlohmann::json>>{krausOps};
      element["probabilities"] = std::vector<double>{1.0};
797
798
      noiseElements.push_back(element);
    }
799
800
801
802
803
  }

  noiseModel["errors"] = noiseElements;
  return noiseModel.dump(6);
}
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835

std::vector<double> IbmqNoiseModel::averageSingleQubitGateFidelity() const {
  std::vector<double> result;
  // Use U3 gate fidelity:
  for (size_t qId = 0; qId < m_nbQubits; ++qId) {
    const std::string gateName = "u3_" + std::to_string(qId);
    const auto gateErrorIter = m_gateErrors.find(gateName);
    const double gateError =
        (gateErrorIter == m_gateErrors.end()) ? 0.0 : gateErrorIter->second;
    result.emplace_back(1.0 - gateError);
  }
  return result;
}

std::vector<std::tuple<size_t, size_t, double>>
IbmqNoiseModel::averageTwoQubitGateFidelity() const {
  std::vector<std::tuple<size_t, size_t, double>> result;
  for (const auto &[gateName, gateError] : m_gateErrors) {
    if (gateName.rfind("cx", 0) == 0) {
      // CNOT gate:
      const std::size_t pos = gateName.find("_");
      const std::string firstArg = gateName.substr(2, pos - 2);
      const std::string secondArg = gateName.substr(pos + 1);
      const auto firstBit = std::atoi(firstArg.c_str());
      const auto secondBit = std::atoi(secondArg.c_str());
      result.emplace_back(
          std::make_tuple(firstBit, secondBit, 1.0 - gateError));
    }
  }

  return result;
}
836
837
838
839
840
841
842
843

std::string
AerAccelerator::getNativeCode(std::shared_ptr<CompositeInstruction> program,
                              const HeterogeneousMap &config) {
  auto qobj_str = xacc_to_qobj->translate(program);
  nlohmann::json j = nlohmann::json::parse(qobj_str)["qObject"];
  return j.dump(2);
}
Mccaskey, Alex's avatar
Mccaskey, Alex committed
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
} // namespace quantum
} // namespace xacc

using namespace cppmicroservices;

namespace {

/**
 */
class US_ABI_LOCAL AerActivator : public BundleActivator {

public:
  AerActivator() {}

  /**
   */
  void Start(BundleContext context) {
    auto xt = std::make_shared<xacc::quantum::AerAccelerator>();
    context.RegisterService<xacc::Accelerator>(xt);
863
864
    context.RegisterService<xacc::NoiseModel>(
        std::make_shared<xacc::quantum::IbmqNoiseModel>());
Mccaskey, Alex's avatar
Mccaskey, Alex committed
865
866
867
868
869
870
871
872
873
874
  }

  /**
   */
  void Stop(BundleContext /*context*/) {}
};

} // namespace

CPPMICROSERVICES_EXPORT_BUNDLE_ACTIVATOR(AerActivator)