#!/usr/bin/env python # coding: utf-8 # In[1]: # Carga de los paquetes que se requieren incialmente. import datetime ejecucion = str(datetime.datetime.now()).replace(":", "_") import sys import ast # permite interpretar la cadena de caracteres como un vector válido dentro del script # In[2]: # Obtener los argumentos pasados al script argumentos = sys.argv # Argumento en el índice 0 es el nombre del script # que coincide con el identificador de mapa # se suprime la extensión del archivo (.py) nombremapa = argumentos[0].rsplit(".", 1)[0] # Restantes argumentos recibidos por el script k = float(argumentos[1]) xc = ast.literal_eval(argumentos[2]) yc = ast.literal_eval(argumentos[3]) token = argumentos[4] backend = argumentos[5] # En ambiente de desarrollo/testing #nombremapa = "mapa00000" # <----- incrementar en cada ejecución #k = 1 #xc = [ 0.51003914, 2.55963937, -0.64283509] #yc = [-4.57068577, 2.98724481, -2.72584275] #token = "ddbcd41b48f7281b2d0acaa742f91055c9b0bf17394ec07922d7964475ff140471a812bc27158e03472edbef5f9f221d452146efdea69f525117347bdd8b7785" #backend = "ibm_nairobi" archivolog = nombremapa + "-" + ejecucion archivores = nombremapa + "+" + ejecucion archivoerr = nombremapa + "a" + ejecucion creacion = 'w' agregado = 'a' # In[3]: def registro(modo, archivo, contenido): contenido = contenido + '\n' with open(archivo + ".txt", modo) as file: file.write(contenido) registro(agregado, archivolog, "xc = " + ', '.join(str(elemento) for elemento in xc) + '\n') registro(agregado, archivolog, "yc = " + ', '.join(str(elemento) for elemento in yc) + '\n') # In[4]: if sys.version_info < (3, 6): registro(creacion, archivoerr, "Utilice Python versión 3.6 o superior.") raise Exception("Utilice Python versión 3.6 o superior.") # In[5]: registro(creacion, archivolog, "Fecha y hora comienzo de la ejecución:" + str(datetime.datetime.now()) + '\n') # In[6]: import numpy as np import matplotlib.pyplot as plt import operator # In[7]: # Paquetes Qiskit from qiskit import IBMQ, assemble from qiskit.providers.ibmq import least_busy from qiskit_optimization import QuadraticProgram from qiskit_optimization.algorithms import MinimumEigenOptimizer from qiskit import BasicAer from qiskit.quantum_info import Pauli from qiskit.utils import QuantumInstance, algorithm_globals from qiskit.algorithms import NumPyMinimumEigensolver, VQE from qiskit.circuit.library import TwoLocal from qiskit.algorithms.optimizers import SPSA # In[8]: registro(agregado, archivolog, "Fecha y hora fin carga módulos:" + str(datetime.datetime.now()) + '\n') registro(agregado, archivolog, "IdMapa = " + nombremapa + '\n') registro(agregado, archivolog, "k = " + str(k) + '\n') registro(agregado, archivolog, "n = " + str(len(xc)) + '\n') registro(agregado, archivolog, "xc = " + ' ' + ', '.join(str(elemento) for elemento in xc) + '\n') registro(agregado, archivolog, "yc = " + ' ' + ', '.join(str(elemento) for elemento in yc) + '\n') registro(agregado, archivolog, "token = " + token + '\n') registro(agregado, archivolog, "backend = " + backend + '\n') # In[9]: # Generador ubicaciónes de los nodos class Initializer: def __init__(self, xc, yc): self.xc = xc self.yc = yc def generate_instance(self): xc = self.xc xy = self.yc n = len(xc) instance = np.zeros([n, n]) for ii in range(0, n): for jj in range(ii + 1, n): instance[ii, jj] = (xc[ii] - xc[jj]) ** 2 + (yc[ii] - yc[jj]) ** 2 instance[jj, ii] = instance[ii, jj] return n, instance # In[10]: # Inicialización del problema con las coordenadas a procesar initializer = Initializer(xc, yc) n, instance = initializer.generate_instance() registro(agregado, archivolog, "instance= \n" + str(instance) + '\n') # In[11]: class QuantumOptimizer: def __init__(self, k, n, instance, token, backend): self.K = k self.n = n self.instance = instance self.token = token self.backend = backend def binary_representation(self, x_sol=0): instance = self.instance n = self.n K = self.K A = np.max(instance) * 100 # Un parámetro de la función de costo # en A se tiene el elemento máximo de # la matríz multiplicado por 100 registro(agregado, archivolog, "A = " + str(A) + '\n') # Determine los pesos w # convierte la matriz bidimensional instance en un vector unidimensional # instance_vec que tiene una longitud de n**2 instance_vec = instance.reshape(n**2) registro(agregado, archivolog, "instance_vec= " + str(instance_vec) + '\n') # Se crea una lista w_list a partir de la lista instance_vec. En la lista # w_list se incluyen aquellos elementos de instance_vec que son mayores a cero w_list = [instance_vec[x] for x in range(n**2) if instance_vec[x] > 0] registro(agregado, archivolog, "w_list = " + str(w_list) + '\n') # Se crea un arreglo de numpy de tamaño n * (n - 1) y llena todas sus entradas con ceros. w = np.zeros(n * (n - 1)) registro(agregado, archivolog, "w = " + str(w) + '\n') # Se copian los elementos de w_list a w en el mismo órden for ii in range(len(w_list)): w[ii] = w_list[ii] # Algunas variables que se usarán # se crea la matriz identidad n x n Id_n = np.eye(n) registro(agregado, archivolog, "Id_n = \n" + str(Id_n) + '\n') # Se crea la matriz unidad (n - 1) x (n - 1) Im_n_1 = np.ones([n - 1, n - 1]) registro(agregado, archivolog, "Im_n_1 = \n" + str(Im_n_1) + '\n') # Se crea un vector de longitud n con todos sus elementos iguales a 1 Iv_n_1 = np.ones(n) registro(agregado, archivolog, "Iv_n_1 = " + str(Iv_n_1) + '\n') # Se sustituye el primer elemento del vecto Iv_n_1 con un 0 Iv_n_1[0] = 0 registro(agregado, archivolog, "Iv_n_1 = " + str(Iv_n_1) + '\n') # Se crea un vector de longitud n - 1 con todos los elementos iguales a 1 Iv_n = np.ones(n - 1) registro(agregado, archivolog, "Iv_n_1 = " + str(Iv_n_1) + '\n') # ? Se crea un vector de longitud n con todos los elementos iguales a 0 # ? excepto el primer elemento que vale 1 neg_Iv_n_1 = np.ones(n) - Iv_n_1 registro(agregado, archivolog, "neg_Iv_n_1= " + str(neg_Iv_n_1) + '\n') # se crea una matriz de n filas por n x (n-1) con entradas 0 y 1 en cierta # disposición particular # Esta matriz se utiliza para transformar un vector de longitud n(n-1) en # una matriz de distancias simétrica n x n. v = np.zeros([n, n * (n - 1)]) registro(agregado, archivolog, "v = \n" + str(v) + '\n') for ii in range(n): count = ii - 1 for jj in range(n * (n - 1)): if jj // (n - 1) == ii: count = ii if jj // (n - 1) != ii and jj % (n - 1) == count: v[ii][jj] = 1.0 # Se crea un vector de largo n x (n - 1) con 0 y 1 en una distribución particular # vn representa la suma de las aristas entrantes de cada uno de los nodos del grafo, # en el contexto del problema para el que se está construyendo la matriz v vn = np.sum(v[1:], axis=0) registro(agregado, archivolog, "vn = " + str(vn) + '\n') # Q define las interacciones entre variables # np.kron(Id_n, Im_n_1) calcula el producto tensorial entre # las matrices de identidad de nxn y la matriz de unos de (n - 1) x (n - 1) # np.dot(v.T, v) realiza el producto punto entre la transpuesta de la matriz v # y la matriz v misma. El resultado es una matriz que contiene productos # internos entre los vectores fila de v. # Por último Q es una matriz que contiene la suma de las matrices # resultantes en las 2 operaciones anteriores multiplicada por un esclar A Q = A * (np.kron(Id_n, Im_n_1) + np.dot(v.T, v)) # g define la contribución de las variables individuales g = ( w - 2 * A * (np.kron(Iv_n_1, Iv_n) + vn.T) - 2 * A * K * (np.kron(neg_Iv_n_1, Iv_n) + v[0].T) ) # c es el desplazamiento constante c = 2 * A * (n - 1) + 2 * A * (K**2) try: max(x_sol) # Evalúa la distancia de costo de una representación binaria de una ruta fun = ( lambda x: np.dot(np.around(x), np.dot(Q, np.around(x))) + np.dot(g, np.around(x)) + c ) cost = fun(x_sol) except: cost = 0 return Q, g, c, cost def construct_problem(self, Q, g, c) -> QuadraticProgram: qp = QuadraticProgram() for i in range(n * (n - 1)): qp.binary_var(str(i)) qp.objective.quadratic = Q qp.objective.linear = g qp.objective.constant = c return qp def solve_problem(self, qp): IBMQ.save_account(self.token, overwrite=True) provider = IBMQ.load_account() backend = provider.get_backend(self.backend) try: algorithm_globals.random_seed = 10598 quantum_instance = QuantumInstance( backend, seed_simulator=algorithm_globals.random_seed, seed_transpiler=algorithm_globals.random_seed, ) vqe = VQE(quantum_instance=quantum_instance) optimizer = MinimumEigenOptimizer(min_eigen_solver=vqe) registro(agregado, archivolog, "Fecha y hora envío a IBM Quantum: " + str(datetime.datetime.now()) + '\n') result = optimizer.solve(qp) registro(agregado, archivolog, "Fecha y hora retorno desde IBM Quantum: " + str(datetime.datetime.now()) + '\n') # Se calcula el costo del resultado obtenido _, _, _, level = self.binary_representation(x_sol=result.x) registro(creacion, archivores, "Se encontró nivel mínimo: " + str(level) + '\n') registro(agregado, archivores, "Ruta optima encontrada : " + str(result.x) + '\n') except Exception as e: registro(creacion, archivoerr, "Ha fayado la ejecución del experimento" + str(e) + '\n') print(f"Job failed: {e}") # In[12]: # Cree una instancia de la clase del optimizador cuántico con parámetros: quantum_optimizer = QuantumOptimizer(k, n, instance, token, backend) # In[13]: # Se crea la reresentación binaria del problema try: Q, g, c, binary_cost = quantum_optimizer.binary_representation() registro(agregado, archivolog, "Q = \n" + str(Q) + '\n') registro(agregado, archivolog, "g = " + str(g) + '\n') registro(agregado, archivolog, "c = " + str(c) + '\n') print("Costo binario:", binary_cost) except NameError as e: print("Advertencia: No se ha ejecutado el código anterior necesario.") registro(creacion, archivoerr, "Advertencia: No se ha ejecutado el código anterior necesario. (" + str(e) + ")" + '\n') print(e) # In[14]: # Se crea el problema a resolver qp = quantum_optimizer.construct_problem(Q, g, c) # In[15]: # se se envía a resuelve el problema a la plataforma de IBM Quantum quantum_optimizer.solve_problem(qp) registro(agregado, archivolog, "Fecha y hora finalización experimento: " + str(datetime.datetime.now()) + '\n')