diff --git a/.pylintrc b/.pylintrc index 7e2e1dc7..c9180322 100644 --- a/.pylintrc +++ b/.pylintrc @@ -268,8 +268,7 @@ variable-naming-style=snake_case # Regular expression matching correct variable names. Overrides variable- # naming-style. If left empty, variable names will be checked with the set # naming style. -variable-rgx=[a-z_][a-z0-9_]{1,30}$ - +variable-rgx=^[a-zA-Z_][a-zA-Z0-9_]{1,30}$ [CLASSES] diff --git a/src/qiskit_qec/codes/classic/fibonacci_code.py b/src/qiskit_qec/codes/classic/fibonacci_code.py new file mode 100644 index 00000000..fc0e6256 --- /dev/null +++ b/src/qiskit_qec/codes/classic/fibonacci_code.py @@ -0,0 +1,89 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020 +# +# This code is licensed under the Apache num_logical_bitsicense, Version 2.0. You may +# obtain a copy of this license in the num_logical_bitsICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/num_logical_bitsICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +"""Defines the Fibonacci Code""" +from typing import Optional + +import numpy as np + +from qiskit_qec.exceptions import QiskitQECError + + +class ClassicFibonacciCode: + """Classic Fib Code""" + + def __init__(self, num_logical_bits: int, codeword_init_arr: Optional[np.array] = None): + """Initialize data for Fibonacci Code. + See arXiv:2002.11738 for more details + We define the Fibonacci code on a two-dimensional lattice + such that it encodes k = L logical bits + using n = (L^2)/2 physical bits. + The exact distance of the code, + d > L, is unknown + + Args: + num_logical_bits (int): Size of codeword, must be a power of 2 + codeword_init_arr (Optional[np.array]): Optional array that + serves as the foundation for building up the rest of the codeword + using cellular automaton update rules, as described in + section II, formula 2 of arXiv:2002.11738. Defaults to None. + + Raises: + QiskitQECError: Invalid input value for num_logical_bits + """ + if num_logical_bits < 4 or not self._is_power_of_2(num_logical_bits): + raise QiskitQECError("num_logical_bits must be >= 4 and a power of 2") + self.num_logical_bits = num_logical_bits + self.codeword_init_arr = self._gen_codeword_init_arr(codeword_init_arr) + self.code_word = self._generate_fibonacci_code_word() + + def _is_power_of_2(self, n): + return n > 0 and (n & (n - 1)) == 0 + + def _gen_codeword_init_arr(self, codeword_init_arr): + if codeword_init_arr is None: + codeword_init_arr = np.zeros(self.num_logical_bits, dtype=int) + codeword_init_arr[((self.num_logical_bits - 1) // 2)] = 1 + return codeword_init_arr + + def _generate_fibonacci_code_word(self): + """Fibonacci Code Codeword Generator. + Generates a valid configuration of classical data bits that can be decoded + by ClassicFibonacciSpanningErrorDecoder and can undergo spanning errors (generate_spanning_error) + + See arXiv:2002.11738 section II figure 3 for more details + """ + # generates from bottom row up + rect_board = np.zeros((self.num_logical_bits // 2, self.num_logical_bits), dtype=np.int8) + rect_board[(self.num_logical_bits // 2) - 1] = self.codeword_init_arr + for row in range((self.num_logical_bits // 2) - 2, -1, -1): + for bit in range(self.num_logical_bits): + new_val = ( + rect_board[row + 1][(bit - 1) % self.num_logical_bits] + ^ rect_board[row + 1][(bit) % self.num_logical_bits] + ^ rect_board[row + 1][(bit + 1) % self.num_logical_bits] + ) + rect_board[row][bit] = new_val + return rect_board + + def __str__(self) -> str: + """Formatted string.""" + return ( + f"Fibonacci Code:\nnum_logical_bits:{self.num_logical_bits}" + f"\ninit_code_word_array: {self.codeword_init_arr}" + ) + + def __repr__(self) -> str: + """String representation.""" + return ( + f"Fibonacci Code:\nnum_logical_bits:{self.num_logical_bits}\ninit_code_word_array:" + f"\n{self.codeword_init_arr}\ncodeword:\n{self.code_word}" + ) diff --git a/src/qiskit_qec/decoders/classic/fibonacci_decoder.py b/src/qiskit_qec/decoders/classic/fibonacci_decoder.py new file mode 100644 index 00000000..fcbde9f4 --- /dev/null +++ b/src/qiskit_qec/decoders/classic/fibonacci_decoder.py @@ -0,0 +1,631 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020 +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +"""Decoder for Spanning Errors on the the ClassicFibonacciCode""" +import copy +import logging +import math +from typing import Dict, Optional, Set, Tuple + +import numpy as np +import pymatching as pm +import rustworkx as rx + +from qiskit_qec.exceptions import QiskitQECError + +logger = logging.getLogger(__name__) + + +class ClassicFibonacciSpanningErrorDecoder: + """ + self.board = + [L*(L//2 - 1).... ((L**2)//2) - 1] + . + . + . + 2L + L + 0 1 2 3 .... L - 1] + """ + + def __init__( + self, + original_errorword: np.array, + halt: int = 9, + ): + """Decoder for Fibonacci Code (ClassicFibonacciCode) + + This decoder is a Minimum Weight Perfect Matching (MWPM) + Decoder implemented in Pymatching, modified to work for + spanning errors on the Fibonacci Code. + See arXiv:2105.13082 for more information about PyMatching. + See section III of arXiv:2002.11738 for more information + about modifying MWPM for spanning errors on the Fibonacci Code. + + L is the number of encoded logical bits + + The code is envisioned as being on a ((L//2, L)) codeboard, + which is represented by an (L//2, L) numpy array. + + self.board is where the results of the decoder + are applied during self.decode_fib_code(). + + Args: + original_errorword (np.array): An (L/bl/2, L) np.array + representing an errored Fibonnaci Code + halt (int, optional): Number of iterations to run + before terminating unsuccessfully. Defaults to 9. + """ + self.num_logical_bits = len(original_errorword[0]) # len + assert ( + math.log2(self.num_logical_bits) % 1 == 0 + ), "L must be some 2**n where n is an int >= 1" + self.no_cols = self.num_logical_bits + self.no_rows = self.num_logical_bits // 2 + self.no_bits = (self.num_logical_bits**2) // 2 # no bits + self.halt = halt + + # fund_sym + self.original_errorword = original_errorword + + self.board = copy.deepcopy(self.original_errorword) + self.board.shape = self.no_bits + self.fund_symmetry = self._generate_init_symmetry() + self.fund_symmetry.shape = (self.no_rows, self.no_cols) + ( + self.fund_stab_parity_check_matrix, + self.fund_parity_rows_to_faces, + ) = self.generate_check_matrix_from_faces(self.fund_symmetry) + self.fund_symmetry.shape = self.no_bits + + self.fund_single_error_syndromes = self.fib_code_fault_enumeration( + self.fund_stab_parity_check_matrix + ) + self.hx_mat = self._generate_plus_x_trans_matrix() + self.hy_mat = self._generate_plus_y_trans_matrix() + + self.fund_stab_faces = copy.deepcopy(self.fund_symmetry) + self.fund_hori_probe_indx = self.no_bits - (self.num_logical_bits // 2) - 1 + self.fund_verti_probe_indx = self.no_bits - 1 + self.fund_stab_faces.shape = ( + self.num_logical_bits // 2, + self.num_logical_bits, + ) # TODO should work + ( + self.fund_check_matrix, + self.board2stab, + ) = self.generate_check_matrix_from_faces(self.fund_stab_faces) + fund_error_pairs = self.fib_code_fault_enumeration(self.fund_check_matrix) + self.fund_matching_graph, self.fundstab2node = self.enumerated_errors2matching_graph( + fund_error_pairs + ) + + self.all_stab_faces = np.ones( + (self.num_logical_bits // 2, self.num_logical_bits), dtype=int + ) + ( + self.all_stabs_check_mat, + self.all_stabs_parity_rows_to_faces, + ) = self.generate_check_matrix_from_faces(self.all_stab_faces) + + # pymatching + self.matching_decoder = pm.Matching(self.fund_matching_graph) + + self.has_decoder_run = False # lets user know whether decoder has already run + + logger.info(" original_errorword is %s", self.original_errorword) + logger.info(" error board is code %s", self.board) + logger.info(" initial symmetry is: %s", self.fund_symmetry) + logger.info( + "fund_stab_parity_check_matrix is: %s", + self.fund_stab_parity_check_matrix, + ) + logger.info("fund_single_error_syndromes is : %s", self.fund_single_error_syndromes) + logger.info(" hx %s", self.hx_mat) + logger.info(" hy is code %s", self.hy_mat) + + def decode_prob(self, syndrome: np.array): + """Returns whether the horizontal probe edge + (aka the bottom middle bit of the probe triangle) + exists an even (0 aka no flip) or odd (1 aka flip) + number of times in the matching graph && the vertical probe + """ + res = self.matching_decoder.decode(syndrome) + return res[self.fund_hori_probe_indx], res[self.fund_verti_probe_indx], res + + def bit_to_rc(self, bit: int): + """Maps maps from bit index (0 to (L**2)//2) - 1) to (row, column) indexing + + In bit notation we think of all the bits of the codeword as being in a line: + So we have bit 0, bit 1, ... all the way until the last bit ((L**2)//2) - 1 + Example: + [0, 1, 2, ................... ((L**2)//2) - 1 ] + + however, we can picture these as being on the + L//2 by L board. + + In numpy, if we reshape the ((L**2)//2) - 1 to a L//2 by L array, + then bit 0 will get mapped to the + (0th row, 0th column). + This is called row, column notation + Bit 1 will get mapped to the 0th row, 1st column. + ... + L//2 gets mapped to the 1st row, 0th column + etc. + + I show here: + + [ [0 1 2 3 ............................... (L//2) -1] + L//2 + . + . + . + [(L - 1) * (L//2) .................... ((L**2)//2) - 1 ] + ] + + Args: + bit (int): Index of the bit in bit notation + """ + row_len = self.num_logical_bits + + rindx = bit // row_len + cindx = bit % row_len + return (rindx, cindx) + + def rc_to_bit(self, row: int, col: int): + """Maps maps from (row, column) index to bit index (0 to (L**2)//2) - 1) + + In bit notation we think of all the bits of the codeword as being in a line: + So we have bit 0, bit 1, ... all the way until the last bit ((L**2)//2) - 1 + Example: + [0, 1, 2, ................... ((L**2)//2) - 1 ] + + however, we can picture these as being on the + L//2 by L board. + + In numpy, if we reshape the ((L**2)//2) - 1 to a L//2 by L array, + then bit 0 will get mapped to the + (0th row, 0th column). + This is called row, column notation + Bit 1 will get mapped to the 0th row, 1st column. + ... + L//2 gets mapped to the 1st row, 0th column + etc. + + I show here: + + [ [0 1 2 3 ............................... (L//2) -1] + L//2 + . + . + . + [(L - 1) * (L//2) .................... ((L**2)//2) - 1 ] + ] + + Args: + row (int): Row index + col (int): Column index + """ + bit = (row * self.num_logical_bits) + col + return bit + + def _generate_plus_x_trans_matrix(self): + """Takes bit to bit + 1 mod rownum aka shifts bit + to the right but wraps around its current row""" + h_mat = np.zeros((self.no_bits, self.no_bits), dtype=int) + for b in range(self.no_bits): + new_bit = self.shift_by_x_scalar(b) + h_mat[new_bit][b] = 1 + return h_mat + + def _generate_plus_y_trans_matrix(self): + """takes bit to (bit + L) mod (L//2) aka shifts bit the row + below but very bottom row shifts to be the 0th row""" + h_mat = np.zeros((self.no_bits, self.no_bits), dtype=int) + for b in range(self.no_bits): + new_bit = self.shift_by_y_scalar(b) + h_mat[new_bit][b] = 1 + return h_mat + + def shift_by_x(self, bitarr: np.array, power: int = 1) -> np.array: + """shifts every entry in array matrix right by power w/ wrap around + + Args: + bitarr (np.array): original code array that serves as a starting point for the shift + Please note, bitarr itself is *not* shifted + power (int, optional): Number of times to shift. Defaults to 1. + + Returns: + np.array: shifted array + """ + + power = power % self.num_logical_bits + hx = np.linalg.matrix_power(self.hx_mat, power) + sol = np.matmul(hx, bitarr) + sol = sol.astype(int) + return sol + + def shift_by_y(self, bitarr: np.array, power: int = 1): + """shifts every entry in array matrix down by power w/ wrap around + + Args: + bitarr (np.array): original code array that serves as a starting point for the shift + Please note, bitarr itself is *not* shifted + power (int, optional): Number of times to shift. Defaults to 1. + + Returns: + np.array: shifted array + """ + + power = power % (self.num_logical_bits // 2) + hy = np.linalg.matrix_power(self.hy_mat, power) + sol = np.matmul(hy, bitarr) + sol = sol.astype(int) + return sol + + def shift_by_y_scalar(self, bit: int, shift_no: int = 1): + """shifts entry in array matrix down by 1 w/ wrap around + + Args: + bit (int): bit to be shifted + shift_no (int, optional): Number of times to shift. Defaults to 1. + + Returns: + int: new location of the bit after shift + """ + + new_bit = bit + for _ in range(shift_no): + new_bit = (new_bit + self.num_logical_bits) % self.no_bits + return new_bit + + def shift_by_x_scalar(self, bit: int, shift_no: int = 1): + """shifts entry in array matrix right by 1 w/ wrap around + Args: + bit (int): bit to be shifted + shift_no (int, optional): Number of times to shift. Defaults to 1. + + Returns: + int: new location of the bit after shift + """ + new_bit = bit + for _ in range(shift_no): + new_bit = ((new_bit + 1) % self.num_logical_bits) + ( + (new_bit // self.num_logical_bits) * (self.num_logical_bits) + ) + return new_bit + + def shift_parity_mat_by_y(self, parity_mat: np.array, power: int = 1): + """Shifts parity matrix vertically by power. + Please note this *modifies* the parity_mat. + + Args: + parity_mat (np.array): parity matrix to be shifted + power (int, optional): Number of times to shift. Defaults to 1. + + Returns: + np.array: original parity_mat + """ + + for row, _ in enumerate(parity_mat): + parity_mat[row] = self.shift_by_y(parity_mat[row], power=power) + return parity_mat + + def shift_parity_mat_by_x(self, parity_mat: np.array, power: int = 1): + """Shifts parity matrix horizontall by power. + Please note this *modifies* the parity_mat. + + Args: + parity_mat (np.array): parity matrix to be shifted + power (int, optional): Number of times to shift. Defaults to 1. + + Returns: + np.array: original parity_mat + """ + for row, _ in enumerate(parity_mat): + parity_mat[row] = self.shift_by_x(parity_mat[row], power=power) + return parity_mat + + def _calc_syndrome(self, check_matr: np.array, board: np.array = None) -> np.array: + """Calculates the syndrome of the check_matrix on the board. + + Args: + check_matr (np.array): check_matrix for calculating the syndrome + board (np.array, optional): board for calculating the syndrome. + If None, defaults to self.board. Defaults to None + + Returns: + np.array: syndrome results + """ + if board is None: + board = self.board + sol = np.matmul(check_matr, board) % 2 + sol = sol.astype(int) + return sol + + def _generate_init_symmetry(self, start_arr: Optional[np.array] = None) -> np.array: + """Generates initial symmetry + + Args: + start_arr (np.array, optional): Optional array that + serves as the foundation for building up the rest of the codeword + using cellular automaton update rules, as described in + section II, formula 2 of arXiv:2002.11738. Defaults to None. + + Raises: + ArgumentError: Validate start_arr + + Returns: + np.array: Initial symmetry of the code + """ + if start_arr and sum(start_arr) != 1: + raise QiskitQECError( + "Can only have a single 1 in start_arr." + + f"All else should be 0 but you have: {start_arr}" + ) + # fund symmetries start from the top instead of the bottom because numpy + rect_board = np.zeros((self.num_logical_bits // 2, self.num_logical_bits), dtype=int) + if start_arr is None: + start_arr = np.zeros(self.num_logical_bits, dtype=int) + start_arr[(self.num_logical_bits // 2) - 1] = 1 + rect_board[0] = start_arr + for row in range(1, self.num_logical_bits // 2): + for bit in range(self.num_logical_bits): + new_val = ( + rect_board[row - 1][(bit - 1) % self.num_logical_bits] + ^ rect_board[row - 1][(bit) % self.num_logical_bits] + ^ rect_board[row - 1][(bit + 1) % self.num_logical_bits] + ) + rect_board[row][bit] = new_val + return rect_board + + def generate_check_matrix_from_faces(self, stab_faces: np.array): + """Use stabilizer faces to make check matrix + + Args: + stab_faces (np.array): Stabilizer faces + + Returns: + np.array: parity matrix + Dict[int, int]: mapping from + parity_check_faces to bit notation, + A dictionary, {x:y}, meaning + that row x of the parity check matrix + has face centered at bit b + + """ + + stab_row2bit_face = {} + parity_mat = np.empty((0, self.no_bits), int) + + stab_row = 0 + for row in range(self.no_rows): + for col in range(self.no_cols): + if stab_faces[row][col] == 1: + a = self.rc_to_bit(row, col) + b = self.rc_to_bit((row) % self.no_rows, (col - 1) % self.no_cols) + c = self.rc_to_bit((row) % self.no_rows, (col + 1) % self.no_cols) + d = self.rc_to_bit( + (row - 1) % self.no_rows, col % self.no_cols + ) # changed to point the other direction + new_stab = np.array([0] * self.no_bits) + new_stab[a] = 1 + new_stab[b] = 1 + new_stab[c] = 1 + new_stab[d] = 1 + parity_mat = np.append(parity_mat, [new_stab], axis=0) + stab_row2bit_face[stab_row] = b + stab_row += 1 + + return parity_mat, stab_row2bit_face + + def fib_code_fault_enumeration( + self, parity_check_matrix: np.array, no_bits: int = None + ) -> Set[Tuple]: + """Return a set containing every possible single + bit flip error and what stabilizers they light up + + Args: + parity_check_matrix (np.array): parity check matrix, + no_bits (int): Number of bits in code array + + Returns: + Set[Tuple]: A set containing every possible single + bit flip error and what stabilizers they light up. + Of the form: {(stab_face1, stab_face2, errorbit0), + (stab_face1, stab_face2, errorbit1), ...} + + """ + if no_bits is None: + no_bits = self.no_bits + error_pairs = set() # (stab_face1, stab_face2, errorbit) + single_error = np.zeros(no_bits, dtype=int) + + for b in range(no_bits): + if no_bits > 10 and b % (no_bits // 10) == 0: + logger.debug("on bit: %s and error set looks like: %s", b, error_pairs) + + ## set up new single error + # clear prev_bit + prev_bit = (b - 1) % no_bits + single_error[prev_bit] = 0 + # set new error + single_error[b] = 1 + + ## what do it light? + lighted = self._calc_syndrome(parity_check_matrix, single_error) + stabs = (lighted == 1).nonzero()[0] + + if len(stabs) % 2 != 0: + emsg = ( + f"Minor panic. Error on bit {b}" + + f"causes a BAD syndrome: {stabs} for lighted: {lighted}" + ) + logger.error( + emsg + ) # TODO just do this via inspection on 1s per column in stab parity check matrix + raise Exception(emsg) + + if len(stabs) > 0: + for indx in range(0, len(stabs), 2): + error_pairs.add((stabs[indx], stabs[indx + 1], b)) + + return error_pairs + + def enumerated_errors2matching_graph( + self, error_graphs: Set[Tuple] + ) -> (rx.PyGraph, Dict[int, int]): + """Turn stabilizer pairs lit by a single error + into edges in the matching graph + + Args: + error_graphs (Set[Tuple]): A set containing + every possible single bit flip error and + what stabilizers they light up. + Of the form: {(stab_face1, stab_face2, errorbit0), + (stab_face1, stab_face2, errorbit1), ...} + + Returns: + rx.PyGraph: Matching Graph for PyMatching + Dict[int, int]: Mapping from stabilizer id to graph node id + + """ + stab2node = {} + graph = rx.PyGraph() + + def add_to_graph(stabid): + if stabid not in stab2node: + nodeid = graph.add_node({"element": stabid}) + stab2node[stabid] = nodeid + return stab2node[stabid] + + for stab0, stab1, fund_e in error_graphs: + graph_node_n0 = add_to_graph(stab0) + graph_node_n1 = add_to_graph(stab1) + graph.add_edge(graph_node_n0, graph_node_n1, {"fault_ids": {fund_e}}) + + return graph, stab2node + + def decode_fib_code(self, force: Optional[bool] = False): + """Run decoder on self.board + + Args: + force (bool, optional): If true, will run + the decoder, even if it has already been + run on self.board + Defaults to False. + + Raises: + QiskitQECError: Raises an exception if + decoder has already been run and force=False + """ + if not force and self.has_decoder_run: + raise QiskitQECError("In order to rerun decoder, you must pass in force=True") + + self.has_decoder_run = True + # generate graphs and mappings + h_correction = np.zeros(self.no_bits, dtype=int) + v_correction = np.zeros(self.no_bits, dtype=int) + parity_check_matrix = copy.deepcopy(self.fund_check_matrix) + hori_probe_indx = self.fund_hori_probe_indx + verti_probe_indx = self.fund_verti_probe_indx + + cur_all_syndrome = prev_all_syndrome = ( + (self._calc_syndrome(self.all_stabs_check_mat, self.board) == 1).sum(), + ) + start_flag = True + meta_round_count = 0 + round_count = 0 + self.fund_stab_faces.shape = self.no_bits + while ( + (cur_all_syndrome < prev_all_syndrome or start_flag) + and cur_all_syndrome != 0 + and meta_round_count < self.halt + ): + start_flag = False + prev_all_syndrome = cur_all_syndrome + + for _ in range(self.num_logical_bits // 2): # will wrap around to all bits + parity_check_matrix = self.shift_parity_mat_by_y(parity_check_matrix) + self.fund_stab_faces = self.shift_by_y(self.fund_stab_faces) + hori_probe_indx = self.shift_by_y_scalar(hori_probe_indx) + verti_probe_indx = self.shift_by_y_scalar(verti_probe_indx) + + for _ in range(self.num_logical_bits): + parity_check_matrix = self.shift_parity_mat_by_x(parity_check_matrix) + self.fund_stab_faces = self.shift_by_x(self.fund_stab_faces) + hori_probe_indx = self.shift_by_x_scalar(hori_probe_indx) + verti_probe_indx = self.shift_by_x_scalar(verti_probe_indx) + + self.fund_stab_faces.shape = ( + self.num_logical_bits // 2, + self.num_logical_bits, + ) + self.fund_stab_faces.shape = self.no_bits + + cur_syndrome = self._calc_syndrome(parity_check_matrix) + # convert syndrome to node + cur_node_syndrome = [0] * len(cur_syndrome) + for stabindx, value in enumerate(cur_syndrome): + nodeindx = self.fundstab2node[stabindx] + cur_node_syndrome[nodeindx] = value # TODO is right? + hcorval, vcorval, res = self.decode_prob(cur_node_syndrome) + + h_correction[hori_probe_indx] = hcorval + v_correction[verti_probe_indx] = vcorval + + round_count += 1 + + logger.debug("current fund stabilizer faces:\n%s", self.fund_stab_faces) + logger.debug("current_parity_check_mat:\n%s", parity_check_matrix) + logger.debug("cur-syndrome-symm: %s", cur_syndrome) + logger.debug("res is: %s", res) + logger.debug("hcorval: %s\nvcorval:%s", hcorval, vcorval) + logger.debug( + "hori probd inex: %s, verti_probe_inx: %s", + hori_probe_indx, + verti_probe_indx, + ) + logger.debug("h_corr: %s\nv_corr:%s", h_correction, v_correction) + + logger.info("Meta-Round %s:", meta_round_count) + logger.info("h_correction: %s\nv_correction:%s", h_correction, v_correction) + + meta_round_count += 1 + d_correction = h_correction * v_correction + hboard = self.board ^ h_correction # apply correction + vboard = self.board ^ v_correction + dboard = self.board ^ d_correction + + hcorsynd = [ + (self._calc_syndrome(self.all_stabs_check_mat, hboard) == 1).sum(), + hboard, + "hori", + ] + vcorsynd = [ + (self._calc_syndrome(self.all_stabs_check_mat, vboard) == 1).sum(), + vboard, + "verti", + ] + dcorsynd = [ + (self._calc_syndrome(self.all_stabs_check_mat, dboard) == 1).sum(), + dboard, + "dcor", + ] + + opts = [hcorsynd, vcorsynd, dcorsynd] + + winner = min(opts, key=lambda x: x[0]) + cur_all_syndrome = winner[0] + self.board = winner[1] # update board to best one + logger.info("Updated board is: \n%s", self.board) + + logger.info("decoding finished") diff --git a/src/qiskit_qec/noise/classic/fibonacci_spanning_noise_generator.py b/src/qiskit_qec/noise/classic/fibonacci_spanning_noise_generator.py new file mode 100644 index 00000000..8ef1e82f --- /dev/null +++ b/src/qiskit_qec/noise/classic/fibonacci_spanning_noise_generator.py @@ -0,0 +1,89 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020 +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +"""Spanning Error Generator for Fibonacci Code""" +import logging +from typing import Optional + +import numpy as np + +from qiskit_qec.codes.classic.fibonacci_code import ClassicFibonacciCode +from qiskit_qec.exceptions import QiskitQECError + +logger = logging.getLogger(__name__) + + +def generate_spanning_error( + code: ClassicFibonacciCode, + offset: int = 0, + probability_of_error: Optional[float] = 1e-7, + width: Optional[int] = None, + is_vertical: Optional[bool] = True, +) -> (np.array, np.array): + """Generates spanning errors for fibonacci code (ClassicFibonacciCode) + that can be fed to the fibonacci decoder (ClassicFibonacciSpanningErrorDecoder) + See section III subsection D figure 7 on arXiv:2002.11738 for more details + + Args: + code (ClassicFibonacciCode): Fibonacci Code, has information about the codeword that + generate_spanning_error will use as a base for it's error + offset (int, optional): Offset for where spanning error begins. Offset will be a + horizontal offset from the left if is_vertical=False and, otherwise, + a vertical offsetfrom the top if is_vertical=True). Defaults to 0. + probability_of_error (float, optional): Error probability. Defaults to 1e-7. + width (int, optional): How wide the error should be. If is_vertical=True, width + willbe how wide the spanning error is across from left to right. Otherwrise, if + is_vertical=False width will be the height of the spanning error from top to bottom. + Defaults to None. + is_vertical (bool, optional): If true, the error will span top to bottom. Otherwise, + the error will span left to right. Defaults to True. + + Raises: + QiskitQECError: probability_of_error must be <= 1 and >=0 + Returns: + error_board (np.array): element-wise "bitwise XOR" of error_mask and codeword + error_mask (np.array): an np.array indicating which bits are flipped by errors + """ + + codeword = code.code_word + if not width: + width = code.num_logical_bits + if probability_of_error < 0 or probability_of_error > 1: + raise QiskitQECError( + f"probability_of_error must be > 0 and < 1 but is: {probability_of_error}" + ) + error_mask = np.zeros(codeword.shape, dtype=int) + + num_rows = len(codeword) + num_col = len(codeword[0]) + if is_vertical: + if probability_of_error == 1: + error_mask[:, offset : width + offset] = 1 + else: + for i in range(num_rows): + for j in range(width): + error_mask[i][(j + offset) % num_col] = np.random.choice( + [0, 1], p=[1 - probability_of_error, probability_of_error] + ) + + else: + height = width + if probability_of_error == 1: + error_mask[offset : offset + height, :] = 1 + else: + for i in range(width): + for j in range(num_col): + error_mask[(i + offset) % num_rows][j] = np.random.choice( + [0, 1], p=[1 - probability_of_error, probability_of_error] + ) + + error_board = error_mask ^ codeword + return error_board, error_mask diff --git a/test/classic/test_fibonacci_code.py b/test/classic/test_fibonacci_code.py new file mode 100644 index 00000000..f2c8f987 --- /dev/null +++ b/test/classic/test_fibonacci_code.py @@ -0,0 +1,57 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020 +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +import unittest + +from qiskit_qec.codes.classic.fibonacci_code import ClassicFibonacciCode +from qiskit_qec.decoders.classic.fibonacci_decoder import ( + ClassicFibonacciSpanningErrorDecoder, +) +from qiskit_qec.noise.classic.fibonacci_spanning_noise_generator import ( + generate_spanning_error, +) + + +class TestClassicFibCode(unittest.TestCase): + """Test Classic Fib Code""" + + def test_full_fib_decoding(self): + test_cases = [ + [4, [(0.01, 1.0), (0.105, 0.8), (0.2, 0.51)]], + [8, [(0.01, 1.0), (0.105, 0.84), (0.2, 0.43)]], + ] + num_shots = 100 + + for test in test_cases: + L = test[0] + for p, accuracy_sol in test[1]: + success_no = 0 + for round in range(num_shots): + code = ClassicFibonacciCode( + L + ) # generate an initial codeword. The default one bit at bottom center and cellular automata rules upward + + error_board, error_mask = generate_spanning_error( + code, probability_of_error=p + ) # setting width to L and vertical=True makes iid noise + f = ClassicFibonacciSpanningErrorDecoder( + error_board + ) # give this class the errored codeword + f.decode_fib_code() + f.board.shape = (L // 2, L) + if (f.board == code.code_word).all(): + success_no += 1 + p_success = success_no / num_shots + self.assertAlmostEqual( + p_success, + accuracy_sol, + delta=0.1, # small decoders are pretty vulnerable to variations in success rates + ), f"L={L} with p={p} gave: {p_success} instead of expected: {accuracy_sol}"