QppAccelerator.cpp 17.3 KB
Newer Older
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh 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 "QppAccelerator.hpp"
14
15
#include <thread>
#include <mutex>
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
16
#include "IRUtils.hpp"
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

namespace {
    inline bool isMeasureGate(const xacc::InstPtr& in_instr)
    {
        return (in_instr->name() == "Measure");
    }

    inline size_t getNumberOfThreads()
    {
        static const size_t NB_THREADS = std::thread::hardware_concurrency();
        return NB_THREADS;
    }

    inline double generateRandomProbability() 
    {
        auto randFunc = std::bind(std::uniform_real_distribution<double>(0, 1), std::mt19937(std::chrono::high_resolution_clock::now().time_since_epoch().count()));
        return randFunc();
    }

    bool applyMeasureOp(KetVectorType& io_psi, size_t in_qubitIndex)
    {
        const auto N = io_psi.size();
        const auto k_range = 1ULL << in_qubitIndex;        
        const double randProbPick = generateRandomProbability();
        double cumulativeProb = 0.0;
        size_t stateSelect = 0;
        //select a state based on cumulative distribution
        while (cumulativeProb < randProbPick && stateSelect < N)
        {
            cumulativeProb += std::norm(io_psi[stateSelect++]);
        }
        stateSelect--;

        //take the value of the measured bit
        bool result = ((stateSelect >> in_qubitIndex) & 1);
        // Collapse the state vector according to the measurement result
        double measProb = 0.0;
        for (size_t g = 0; g < N; g += (k_range * 2))
        {
            for (size_t i = 0; i < k_range; ++i)
            {
                if ((((i + g) >> in_qubitIndex) & 1) == result)
                {
                    measProb += std::norm(io_psi[i + g]);
                    io_psi[i + g + k_range] = 0.0; 
                }
                else
                {
                    measProb += std::norm(io_psi[i + g + k_range]);
                    io_psi[i + g] = 0.0;
                }
            }
        }
        // Renormalize the state
        measProb = std::sqrt(measProb);
        for (size_t g = 0; g < N; g += (k_range * 2))
        {
            for (size_t i = 0; i < k_range; i += 1)
            {
                if ((((i + g) >> in_qubitIndex) & 1) == result)
                {
                    io_psi[i + g] /= measProb;
                }
                else
                {
                    io_psi[i + g + k_range] /= measProb;
                }    
            }
        }

        return result;
    }

    void generateMeasureBitString(std::shared_ptr<xacc::AcceleratorBuffer> in_buffer, const std::vector<size_t>& in_bits, const KetVectorType& in_stateVec, int in_shotCount, bool in_multiThread = false)
    {
        if (!in_multiThread)
        {
            // Sequential execution
            for (int i = 0; i < in_shotCount; ++i)
            {
                std::string bitString;
                auto stateVecCopy = in_stateVec;
                for (const auto& bit : in_bits)    
                {
                    bitString.append(std::to_string(applyMeasureOp(stateVecCopy, bit)));
                }

                in_buffer->appendMeasurement(bitString);
            }
        }
        else
        {
            // Parallel execution: divide the shot counts to threads.
            std::vector<std::string> bitStringArray(in_shotCount);
            std::vector<std::thread> threads(getNumberOfThreads());
            std::mutex critical;
            for(int t = 0; t < getNumberOfThreads(); ++t)
            {
                threads[t] = std::thread(std::bind([&](int beginIdx, int endIdx, int threadIdx) {
                    for(int i = beginIdx; i < endIdx; ++i)
                    {
                        std::string bitString;
                        auto stateVecCopy = in_stateVec;
                        for (const auto& bit : in_bits)    
                        {
                            bitString.append(std::to_string(applyMeasureOp(stateVecCopy, bit)));
                        }
                        bitStringArray[i] = bitString;
                    }
                    {
                        // Add measurement bitstring to the buffer:
                        std::lock_guard<std::mutex> lock(critical);
                        for(int i = beginIdx; i < endIdx; ++i)
                        {
                            in_buffer->appendMeasurement(bitStringArray[i]);
                        }
                    }
                }, 
                t * in_shotCount / getNumberOfThreads(), 
                (t+1) == getNumberOfThreads() ? in_shotCount: (t+1) * in_shotCount/getNumberOfThreads(), 
                t));
            }
            std::for_each(threads.begin(),threads.end(),[](std::thread& x){
                x.join();
            });
        }
    }

    // Helper to determine if shot count distribution can be simulated by 
    // sampling the final state vector.
147
    // Requirements: no reset and measure at the very end and nothing after measure.
148
149
150
    bool shotCountFromFinalStateVec(const std::shared_ptr<CompositeInstruction>& in_composite)
    {
        InstructionIterator it(in_composite);
151
        std::set<std::size_t> bitsMeasured;
152
        bool hasReset = false;
153
        bool postMeasureGates = false;
154
155
156
157
158
        while (it.hasNext())
        {
            auto nextInst = it.next();
            if (nextInst->isEnabled())
            {
159
160
                if (nextInst->name() == "Reset") {
                    hasReset = true;
161
                    break;
162
163
                }

164
165
                auto bits = nextInst->bits();

166
167
                if (isMeasureGate(nextInst))
                {
168
                    bitsMeasured.insert(bits.begin(), bits.end());
169
                }
170
                else
171
                {
172
173
174
175
176
177
178
179
                    postMeasureGates = std::any_of(bits.begin(), bits.end(),
                                                   [&](std::size_t bit) {
                                                       return bitsMeasured.count(bit);
                                                   });
                    if (postMeasureGates)
                    {
                        break;
                    }
180
181
182
183
                }
            }
        }

184
        // If Measure gates are at the very end and no reset,
185
        // this Composite can be simulated by random sampling from the state vec.
186
        return !hasReset && !postMeasureGates;
187
    }
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216

    Eigen::MatrixXcd convertToEigenMat(const NoiseModelUtils::cMat& in_stdMat)
    {
        Eigen::MatrixXcd result =  Eigen::MatrixXcd::Zero(in_stdMat.size(), in_stdMat.size());
        for (size_t row = 0; row < in_stdMat.size(); ++row)
        {
            for (size_t col = 0; col < in_stdMat.size(); ++col)
            {
                result(row, col) = in_stdMat[row][col];
            }
        }
        return result;
    }

    NoiseModelUtils::cMat convertToStdMat(const Eigen::MatrixXcd& in_eigenMat)
    {
        const size_t dim = in_eigenMat.rows();
        NoiseModelUtils::cMat result;
        for (size_t row = 0; row < dim; ++row)
        {
            std::vector<std::complex<double>> rowVec;
            for (size_t col = 0; col < dim; ++col)
            {
                rowVec.emplace_back(in_eigenMat(row, col)); 
            }
            result.emplace_back(rowVec);
        }
        return result;
    }
217
}
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
218
219
220
221

namespace xacc {
namespace quantum {

222
    void QppAccelerator::initialize(const HeterogeneousMap& params)
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
223
224
    {
        m_visitor = std::make_shared<QppVisitor>();
225
226
        // Default: no shots (unless otherwise specified)
        m_shots = -1;
227
        if (params.keyExists<int>("shots"))
228
229
        {
            m_shots = params.get<int>("shots");
230
            if (m_shots < 1)
231
232
233
234
            {
                xacc::error("Invalid 'shots' parameter.");
            }
        }
235
236
237
        // Enable VQE mode by default if not using shots.
        // Note: in VQE mode, only expectation values are computed.
        m_vqeMode = (m_shots < 1);
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
238
239
240
241
242
243
244
245
        if (params.keyExists<bool>("vqe-mode"))
        {
            m_vqeMode = params.get<bool>("vqe-mode");
            if (m_vqeMode)
            {
                xacc::info("Enable VQE Mode.");
            }
        }
246
247
248
249

        if (params.keyExists<std::vector<std::pair<int,int>>>("connectivity")) {
            m_connectivity = params.get<std::vector<std::pair<int,int>>>("connectivity");
        }
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
250
251
    }

252
    void QppAccelerator::execute(std::shared_ptr<AcceleratorBuffer> buffer, const std::shared_ptr<CompositeInstruction> compositeInstruction)
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
253
    {
254
255
        const auto runCircuit = [&](bool shotsMode){
            m_visitor->initialize(buffer, shotsMode);
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
256

257
258
            // Walk the IR tree, and visit each node
            InstructionIterator it(compositeInstruction);
259
            while (it.hasNext())
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
260
            {
261
                auto nextInst = it.next();
262
                if (nextInst->isEnabled())
263
                {
264
265
266
267
268
269
                    try {
                        nextInst->accept(m_visitor);
                    } catch (std::exception& ex) {
                        std::cout <<"  QPP CAUGHT EXCEPTION:\n";
                        xacc::error("");
                    }
270
                }
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
271
272
            }

273
274
275
            m_visitor->finalize();
        };

276
277
278
        // Not possible to simulate shot count by direct sampling,
        // e.g. must collapse the state vector.
        if(!shotCountFromFinalStateVec(compositeInstruction))
279
        {
280
281
282
283
284
285
            if (m_shots < 0)
            {
                runCircuit(false);
            }
            else
            {
286
287
288
289
                for (int i = 0; i < m_shots; ++i)
                {
                    runCircuit(true);
                }
290
291
292
293
294
295
296
297
298
299
300
301
302
            }
        }
        else
        {
            // Index of measure bits
            std::vector<size_t> measureBitIdxs;
            m_visitor->initialize(buffer);
            // Walk the IR tree, and visit each node
            InstructionIterator it(compositeInstruction);
            while (it.hasNext())
            {
                auto nextInst = it.next();
                if (nextInst->isEnabled())
303
                {
304
                    if (!isMeasureGate(nextInst))
305
                    {
306
307
308
309
310
311
                        nextInst->accept(m_visitor);
                    }
                    else
                    {
                        // Just collect the indices of measured qubit
                        measureBitIdxs.emplace_back(nextInst->bits()[0]);
312
313
                    }
                }
314
315
316
317
318
319
320
321
322
323
324
325
            }
            
            // Run bit-string simulation
            if (!measureBitIdxs.empty())
            {
                const auto& stateVec = m_visitor->getStateVec();
                if (m_shots < 0)
                {
                    const double expectedValueZ = QppVisitor::calcExpectationValueZ(stateVec, measureBitIdxs);
                    buffer->addExtraInfo("exp-val-z", expectedValueZ);
                }
                else
326
327
328
329
330
                {
                    // Try multi-threaded execution if there are many shots.
                    const bool multiThreadEnabled = (m_shots > 1024);
                    generateMeasureBitString(buffer, measureBitIdxs, stateVec, m_shots, multiThreadEnabled);
                }
331
            }
332
333
            // Note: must save the state-vector before finalizing the visitor.
            cacheExecutionInfo();
334
            m_visitor->finalize();
335
        }
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
336
    }
337
    
338
    void QppAccelerator::execute(std::shared_ptr<AcceleratorBuffer> buffer, const std::vector<std::shared_ptr<CompositeInstruction>> compositeInstructions)
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
339
    {
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
340
        if (!m_vqeMode || compositeInstructions.size() <= 1) 
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
341
        {
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
            for (auto& f : compositeInstructions)
            {
                auto tmpBuffer = std::make_shared<xacc::AcceleratorBuffer>(f->name(), buffer->size());
                execute(tmpBuffer, f);
                buffer->appendChild(f->name(), tmpBuffer);
            }
        }
        else 
        {
            auto kernelDecomposed = ObservedAnsatz::fromObservedComposites(compositeInstructions);
            // Always validate kernel decomposition in DEBUG
            assert(kernelDecomposed.validate(compositeInstructions));
            m_visitor->initialize(buffer);            
            // Base kernel:
            auto baseKernel = kernelDecomposed.getBase();
            // Basis-change + measures
            auto obsCircuits = kernelDecomposed.getObservedSubCircuits();
            // Walk the base IR tree, and visit each node
360
            InstructionIterator it(baseKernel);
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
            while (it.hasNext()) 
            {
                auto nextInst = it.next();
                if (nextInst->isEnabled() && !nextInst->isComposite()) 
                {
                    nextInst->accept(m_visitor);
                }
            }

            // Now we have a wavefunction that represents execution of the ansatz.
            // Run the observable sub-circuits (change of basis + measurements)
            for (int i = 0; i < obsCircuits.size(); ++i) 
            {
                auto tmpBuffer = std::make_shared<xacc::AcceleratorBuffer>(obsCircuits[i]->name(), buffer->size());
                const double e = m_visitor->getExpectationValueZ(obsCircuits[i]);
                tmpBuffer->addExtraInfo("exp-val-z", e);
                buffer->appendChild(obsCircuits[i]->name(), tmpBuffer);
            }
            m_visitor->finalize();
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
380
381
        }
    }
382

383
384
    void QppAccelerator::apply(std::shared_ptr<AcceleratorBuffer> buffer, std::shared_ptr<Instruction> inst) 
    {
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
385
386
        if (!m_visitor->isInitialized()) {
            m_visitor->initialize(buffer);
387
            m_currentBuffer = std::make_pair(buffer.get(), buffer->size());
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
388
        }
389
390
391
392
393
394
395
396
397
398
399
400
401
402

        if (buffer.get() == m_currentBuffer.first &&
            buffer->size() != m_currentBuffer.second) {
            if (buffer->size() > m_currentBuffer.second) {
                const auto nbQubitsToAlloc = buffer->size() - m_currentBuffer.second;
                // std::cout << "Allocate " << nbQubitsToAlloc << " qubits.\n";
                m_visitor->allocateQubits(nbQubitsToAlloc);
                m_currentBuffer.second = buffer->size();
            }
            else {
                xacc::error("Qubits deallocation is not supported.");
            }
        }

403
404
405
406
407
408
409
        if (inst->isComposite() || inst->isAnalog())
        {
            xacc::error("Only gates are allowed.");
        }
        if (inst->name() == "Measure")
        {
            const auto measRes = m_visitor->measure(inst->bits()[0]);
410
            buffer->measure(inst->bits()[0], (measRes ? 1 : 0));
411
412
413
414
415
        }
        else
        {
            auto gateCast = std::dynamic_pointer_cast<xacc::quantum::Gate>(inst);
            assert(gateCast);
416
417
418
419
420
421
422
423
            try {
                m_visitor->applyGate(*gateCast);
            } catch(std::exception& e) {
                // std::cout << "CAUGHT QPP EXCP\n";
                std::stringstream ss;
                ss << "Error executing QPP::apply() on gate: " << gateCast->toString() << "\n";
                xacc::error(ss.str());
            }
424
425
426
        }
    }

427
428
429
430
431
432
433
434
435
436
437
438
    void QppAccelerator::cacheExecutionInfo() {
      // Cache the state-vector:
      // Note: qpp stores wavefunction in Eigen vectors,
      // hence, maps to std::vector.
      auto stateVec = m_visitor->getStateVec();
      ExecutionInfo::WaveFuncType waveFunc(stateVec.data(),
                                           stateVec.data() + stateVec.size());
      m_executionInfo = {
          {ExecutionInfo::WaveFuncKey,
           std::make_shared<ExecutionInfo::WaveFuncType>(std::move(waveFunc))}};
    }

439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
    NoiseModelUtils::cMat DefaultNoiseModelUtils::krausToChoi(const std::vector<NoiseModelUtils::cMat>& in_krausMats) const
    {
        std::vector<Eigen::MatrixXcd> krausMats;
        for (const auto& mat : in_krausMats)
        {
            krausMats.emplace_back(convertToEigenMat(mat));
        }
        return convertToStdMat(qpp::kraus2choi(krausMats));
    }

    std::vector<NoiseModelUtils::cMat> DefaultNoiseModelUtils::choiToKraus(const NoiseModelUtils::cMat& in_choiMat) const
    {
        std::vector<NoiseModelUtils::cMat> resultKraus;
        const auto krausMats = qpp::choi2kraus(convertToEigenMat(in_choiMat));
        for (const auto& mat : krausMats)
        {
            resultKraus.emplace_back(convertToStdMat(mat));
        }
        return resultKraus;
    }
459
460
461
462
463
464
465
466
467
468
469
470
471
472

    NoiseModelUtils::cMat DefaultNoiseModelUtils::combineChannelOps(const std::vector<NoiseModelUtils::cMat> &in_choiMats) const 
    {
        assert(!in_choiMats.empty());
        auto choiSum = convertToEigenMat(in_choiMats[0]);
        for (size_t i = 1; i < in_choiMats.size(); ++i)
        {
            const auto nextOp = convertToEigenMat(in_choiMats[i]);
            choiSum = choiSum + nextOp;
        }
        const double normalized = std::abs((choiSum(0,0) + choiSum(1,1) + choiSum(2,2) + choiSum(3,3)).real()/2.0);
        choiSum = (1/normalized)*choiSum;
        return convertToStdMat(choiSum);
    }
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
473
}}