qaoa_maxcut.cpp 13.5 KB
Newer Older
1
/*******************************************************************************
2
 * Copyright (c) 2021 UT-Battelle, LLC.
3
4
5
6
7
8
9
10
11
12
13
 * 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:
 *   A.M. Santana - initial API and implementation
 *******************************************************************************/

14
#include "qaoa_maxcut.hpp"
15

16
#include "xacc.hpp"
17
#include "Circuit.hpp"
18
#include "xacc_service.hpp"
19
#include "PauliOperator.hpp"
20
#include "xacc_observable.hpp"
21
22
#include "CompositeInstruction.hpp"
#include "AlgorithmGradientStrategy.hpp"
Nguyen, Thien's avatar
Nguyen, Thien committed
23
#include "IRTransformation.hpp"
24
25
#include <cassert>
#include <iomanip>
26
#include <random>
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
147
148
149
150
151
152


// WARM-START FUNCTIONS:
auto random_vector(const double l_range, const double r_range, const std::size_t size) {
    // Generate a random initial parameter set 
    // WARNING: This function multiplies each output in the range
    // by 2*pi
    std::random_device rnd_device;
    std::mt19937 mersenne_engine{rnd_device()}; // Generates random integers
    std::uniform_real_distribution<double> dist{l_range, r_range};
    auto gen = [&dist, &mersenne_engine]() { return dist(mersenne_engine); };
    std::vector<double> vec(size);
    std::generate(vec.begin(), vec.end(), gen);
    for (int i=0; i < vec.size(); ++i){vec.at(i) *= 2 * M_PI;}
    return vec;
}

double cutValue2D(xacc::Graph* const graph, std::vector<double> positions)
{
    double value = 0.0;
    double w, ux, uy, vx, vy;
    for (int i = 0; i < graph->order(); ++i) {
        auto neighbors = graph->getNeighborList(i);
        for (const auto &neighborId : neighbors) {
            w = graph->getEdgeWeight(i,neighborId) ? graph->getEdgeWeight(i,neighborId) : 1.0;
            ux = std::cos(positions[i]);
            uy = std::sin(positions[i]);
            vx = std::cos(positions[neighborId]);
            vy = std::sin(positions[neighborId]);
            value += w * (std::pow(ux-vx,2) + std::pow(uy-vy,2));
        }
    }
    return (value / 4);
}

std::pair<double,double> regularize(double in_phi, double in_theta)
{
    auto theta = std::fmod(in_theta, 2 * M_PI);
    auto phi = in_phi;
    if (theta > M_PI)
    {
        theta = (2*M_PI) - theta;
        phi = std::fmod(in_phi + M_PI, 2 * M_PI);
    }
    return std::make_pair(phi,theta);
}

std::tuple<double,double,double> polartocart(double phi, double theta) 
{
    auto x = std::sin(theta) * std::cos(phi);
    auto y = std::sin(phi) * std::sin(theta);
    auto z = std::cos(theta);
    return std::make_tuple(x,y,z);
}

std::pair<double,double> carttopolar(double x, double y, double z)
{
    double phi = std::atan2(y,x);
    double theta = std::acos(z / std::sqrt(x*x+y*y+z*z));
    return std::make_pair(phi,theta);
}

std::pair<std::vector<double>, std::vector<double>> getSDPSolution(xacc::Graph* const graph)
{
    int n_nodes = graph->order();
    float W = 0;
    int trials = 1;

    float maxDelta = 10;
    std::vector<double> values, positions;
    double cutvalue, dTheta;
    // Setting best value to -inf.
    const int bestValue = std::numeric_limits<int>::min();
    for (int j=0; j < trials; ++j){
        positions = random_vector(0, 1, n_nodes);
        while (values.size() <= 100) {
            for (int i=0; i < n_nodes; ++i){
                cutvalue = cutValue2D(graph, positions);
                dTheta = (std::rand() - 0.5) / maxDelta;
                positions[i] += dTheta;
                if (cutValue2D(graph, positions) < cutvalue){
                    positions[i] -= dTheta;
                }
            }
            values.push_back(cutvalue);
        }
    }
    std::pair<std::vector<double>, std::vector<double>> solution = std::make_pair(positions, values);
    return solution;
}

std::pair<std::vector<double>, std::vector<double>> rotatesolution(std::pair<std::vector<double>, std::vector<double>> solution, int top)
{   
    std::vector<double> positions = std::get<0>(solution);
    std::vector<double> values = std::get<1>(solution);
    auto angles = regularize(positions[top-1], values[top-1]);
    auto phi = std::get<0>(angles);
    auto theta = std::get<1>(angles);
    theta *= -1.;
    std::vector<double> newPositions, newValues;
    std::tuple<double,double,double> xyz;
    std::pair<double,double> polar;
    double newX, newY, newZ;
    int i = 0;
    for (auto position : positions)
    {
        position = std::fmod(position - phi, 2 * M_PI);
        auto newvals = regularize(position, values[i]);
        position = std::get<0>(newvals);
        values[i] = std::get<1>(newvals);
        xyz = polartocart(position, values[i]);
        newX = std::cos(theta) * std::get<0>(xyz) * std::sin(theta) * std::get<2>(xyz);
        newY = std::get<1>(xyz);
        newZ = -1. * std::sin(theta) * std::get<0>(xyz) * std::cos(theta) * std::get<2>(xyz);
        polar = carttopolar(newX,newY,newZ);
        newPositions.push_back(std::get<0>(polar));
        newValues.push_back(std::get<1>(polar));
        i += 1;
    }
    phi = 2 * M_PI * std::rand();
    for (int i=0; i < newPositions.size(); ++i)
    {
        newPositions[i] = std::fmod(newPositions[i] + phi, 2 * M_PI);
    }
    std::pair<std::vector<double>, std::vector<double>> rotated_solution = std::make_pair(newPositions, newValues);
    return rotated_solution;
153
154
} 
// end warm_start functions
155
156
157

namespace xacc {
namespace algorithm {
158
159
    // Construct the MAX-CUT Hamiltonian as an XACC observable
    // based on input graphs edges
Santana, Anthony's avatar
Santana, Anthony committed
160
    std::shared_ptr<Observable> maxcut_qaoa::constructMaxCutHam(xacc::Graph *in_graph) const {
161
162
163
164
165
166
167
168
169
170
171
172
        xacc::quantum::PauliOperator H;
        for (int i = 0; i < in_graph->order(); ++i) {
            auto neighbors = in_graph->getNeighborList(i);
            for (const auto &neighborId : neighbors) {
                if (in_graph->getEdgeWeight(i,neighborId)){
                    auto weight = in_graph->getEdgeWeight(i,neighborId);
                    H += 0.5 * weight * (xacc::quantum::PauliOperator(1.) - xacc::quantum::PauliOperator({{i, "Z"}, {neighborId, "Z"}}));
                } else {
                H += 0.5 * (xacc::quantum::PauliOperator(1.) - xacc::quantum::PauliOperator({{i, "Z"}, {neighborId, "Z"}})); 
                }
            }
        }
Santana, Anthony's avatar
Santana, Anthony committed
173
        return xacc::quantum::getObservable("pauli", H.toString());
174
175
176
    }

    // Take a graph as an input and output a series of Rx/Rz instructions on 
177
    // each qubit for Warm-Start procedure (https://arxiv.org/abs/2010.14021)
178
179
180
    auto warm_start(xacc::Graph* m_graph) {
        auto provider = getIRProvider("quantum");
        auto initial_state = provider->createComposite("initial_state");
181
        auto solution1 = getSDPSolution(m_graph);
182
183
        int sol_size = std::get<0>(solution1).size() + std::get<1>(solution1).size();
        int rand = std::rand() % sol_size;
184
185
186
        auto rotated = rotatesolution(solution1, rand);
        std::vector<double> x_rotations = std::get<1>(rotated);
        std::vector<double> z_rotations = std::get<0>(rotated);
Santana, Anthony's avatar
Santana, Anthony committed
187
        for (size_t i = 0; i < m_graph->order(); ++i) {
188
189
            auto xrot = x_rotations[i];
            auto zrot = z_rotations[i] + (M_PI / 2);
190
191
            initial_state->addInstruction(provider->createInstruction("Rx", { i }, {xrot}));
            initial_state->addInstruction(provider->createInstruction("Rz", { i }, {zrot}));
192
        }
193
        return initial_state;
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
    }
    
    bool maxcut_qaoa::initialize(const HeterogeneousMap &parameters) {
        bool initializeOk = true;
        // Hyper-parameters for QAOA:
        // (1) Accelerator (QPU)
        if (!parameters.pointerLikeExists<Accelerator>("accelerator")) {
            std::cout << "'accelerator' is required.\n";
            // We check all required params; hence don't early return on failure.
            initializeOk = false;
        }

        // (2) Classical optimizer
        if (!parameters.pointerLikeExists<Optimizer>("optimizer")) {
            std::cout << "'optimizer' is required.\n";
            initializeOk = false;
        }

        // (3) Number of mixing and cost function steps to use (default = 1)
        m_nbSteps = 1;
        if (parameters.keyExists<int>("steps")) {
            m_nbSteps = parameters.get<int>("steps");
        }

        // (4) Cost Hamiltonian or a graph to construct the max-cut cost Hamiltonian
        // from.
        bool graphInput = false;
        if (parameters.pointerLikeExists<Graph>("graph")) {
            graphInput = true;
            m_graph = parameters.getPointerLike<Graph>("graph");
        } else {
225
            std::cout << "graph' is required.\n";
226
227
228
229
230
231
232
233
234
235
            initializeOk = false;
        } 
        
        // Default is Extended ParameterizedMode (less steps, more params)
        m_parameterizedMode = "Extended";
        if (parameters.stringExists("parameter-scheme")) {
            m_parameterizedMode = parameters.getString("parameter-scheme");
        }

        if (initializeOk) {
Santana, Anthony's avatar
Santana, Anthony committed
236
            m_costHamObs = constructMaxCutHam(parameters.getPointerLike<Graph>("graph"));
237
238
239
240
241
242
243
            m_qpu = parameters.getPointerLike<Accelerator>("accelerator");
            m_optimizer = parameters.getPointerLike<Optimizer>("optimizer");
            // Optional ref-hamiltonian
            m_refHamObs = nullptr;
            if (parameters.pointerLikeExists<Observable>("ref-ham")) {
            m_refHamObs = parameters.getPointerLike<Observable>("ref-ham");
            }
Nguyen, Thien's avatar
Nguyen, Thien committed
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
            m_irTransformation = nullptr;
            // This QPU has topology-constraint
            if (!m_qpu->getConnectivity().empty()) {
              if (parameters.pointerLikeExists<xacc::IRTransformation>(
                      "placement")) {
                m_irTransformation = xacc::as_shared_ptr(
                    parameters.getPointerLike<xacc::IRTransformation>(
                        "placement"));
                if (m_irTransformation->type() !=
                    xacc::IRTransformationType::Placement) {
                  xacc::error(m_irTransformation->name() +
                              " is not a placement service.");
                }
              }
            }
259
260
261
262
263
264
265
266
267
268
269
        }

        // Check if a parameter initialization routine has been specified
        // Current options: random (default), warm-starts, (eventually) Fourier 
        if (parameters.stringExists("initialization")) {
            m_initializationMode = parameters.getString("initialization");
        }

        if (parameters.keyExists<bool>("maximize")) {
            m_maximize = parameters.get<bool>("maximize");
        }
270
271
272
273
        m_shuffleTerms = false;
        if (parameters.keyExists<bool>("shuffle-terms")) {
          m_shuffleTerms = parameters.get<bool>("shuffle-terms");
        }
274
275
276
277
278
        if (m_optimizer && m_optimizer->isGradientBased() &&
            gradientStrategy == nullptr) {
            // No gradient strategy was provided, just use autodiff.
            gradientStrategy = xacc::getService<AlgorithmGradientStrategy>("autodiff");
            gradientStrategy->initialize(
Santana, Anthony's avatar
Santana, Anthony committed
279
                {{"observable", m_costHamObs}});
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
        }

        return initializeOk;
        }

        const std::vector<std::string> maxcut_qaoa::requiredParameters() const {
        return {"accelerator", "optimizer", "observable"};
        }
    
    void maxcut_qaoa::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
        HeterogeneousMap m; 
        m.insert("accelerator", m_qpu);
        m.insert("optimizer", m_optimizer);
        m.insert("observable", m_costHamObs);
        m.insert("steps", m_nbSteps);
        m.insert("parameter-scheme", "Standard");
296
        m.insert("shuffle-terms", m_shuffleTerms);
297
298

        if (m_initializationMode == "warm-start") {
299
300
            m.insert("initial-state", warm_start(m_graph));
        } else if (m_initial_state){
301
            m.insert("initial-state", m_initial_state);
302
        }
Nguyen, Thien's avatar
Nguyen, Thien committed
303
304
305
        if (m_irTransformation) {
            m.insert("placement", m_irTransformation);
        }
306
307
308
309
310
311
        // Initialize QAOA
        auto qaoa = xacc::getAlgorithm("QAOA", m);
        // Allocate some qubits and execute
        qaoa->execute(buffer);
    }

312
    std::vector<double> maxcut_qaoa::execute(const std::shared_ptr<AcceleratorBuffer> buffer, const std::vector<double> &x) {
313
        // Build up heterogeneous map of algorithm options to pass off to "QAOA" algorithm
314
315
316
317
318
319
        HeterogeneousMap m; 
        m.insert("accelerator", m_qpu);
        m.insert("optimizer", m_optimizer);
        m.insert("observable", m_costHamObs);
        m.insert("steps", m_nbSteps);
        m.insert("parameter-scheme", "Standard");
320
321
        m.insert("shuffle-terms", m_shuffleTerms);

322
323
324
325
        // If: warm-start selected, call the warm start function
        // and use as initial-state.
        // Else if: custom initial-state given, pass that to the
        // algorithm call.
326
        // Else: "QAOA" algorithm defaults to Hadamards on each qubit.
327
328
329
330
331
332
        if (m_initializationMode == "warm-start") {
            m.insert("initial-state", warm_start(m_graph));
        } else if (m_initial_state){
            m.insert("initial-state", m_initial_state);
        }

Nguyen, Thien's avatar
Nguyen, Thien committed
333
334
335
336
        if (m_irTransformation) {
            m.insert("placement", m_irTransformation);
        }

337
338
339
340
341
342
343
        // Initialize QAOA
        auto qaoa = xacc::getAlgorithm("QAOA", m);
        // Allocate some qubits and execute
        auto energy = qaoa->execute(buffer, x);
        return energy;
    }

344
345
} // namespace quantum
} // namespace xacc