from qiskit import ClassicalRegister, QuantumRegister
from qiskit import QuantumCircuit, execute
from qiskit import register, available_backends, get_backend

#import Qconfig and set APIToken and API url
import sys
sys.path.append("../../") # go to parent dir
import Qconfig
qx_config = {
    "APItoken": Qconfig.APItoken,
    "url": Qconfig.config['url']}

#set api
register(qx_config['APItoken'], qx_config['url'])

import math, scipy, copy, time
import numpy as np

def AddDelay (script,q,num,delay):
    
    for _ in range(delay):
        for address in range(num):
            script.iden(q[address])

def  AddCnot(script,q,control,target,device):
    
    # set the coupling map
    # b in coupling_map[a] means a CNOT with control qubit a and target qubit b can be implemented
    # this is essentially just copy and pasted from
    #     https://github.com/IBM/qiskit-qx-info/tree/master/backends/ibmqx3
    #     https://github.com/IBM/qiskit-qx-info/tree/master/backends/ibmqx5
    # but with a few stylistic changes
    coupling_map = {}
    coupling_map['ibmqx3'] = {0: [1], 1: [2], 2: [3], 3: [14], 4: [3, 5], 5: [], 6: [7, 11], 7: [10],
                              8: [7], 9: [10, 8], 10:[], 11: [10], 12: [5, 11, 13], 13: [4, 14], 14:[], 15: [0, 14]}
    coupling_map['ibmqx5'] = {0:[],1:[0,2], 2:[3], 3:[4, 14], 4:[], 5:[4], 6:[5,7,11], 7:[10],
                              8:[7], 9:[8, 10], 10:[], 11:[10], 12:[5, 11, 13], 13:[4, 14], 14:[], 15:[0, 2, 14]}
        
    simulator = (device not in ['ibmqx3','ibmqx5'])
    
    # if simulator, just do the CNOT
    if simulator:
        script.cx(q[control], q[target])
    # same if the coupling map allows it
    elif target in coupling_map[device][control]:
        script.cx(q[control], q[target])
    # if it can be done the other way round we conjugate with Hadamards
    elif ( control in coupling_map[device][target] ):
        script.h(q[control])
        script.h(q[target])
        script.cx(q[target], q[control])
        script.h(q[control])
        script.h(q[target])
    else:
        print('Qubits ' + str(control) + ' and ' + str(target) + ' cannot be entangled.')

def GetAddress (codeQubit,offset,simulator):
    
    if (simulator):
        address = 2*codeQubit + offset
    else:
        address = (5-2*codeQubit-offset)%16
        #address = (6+2*codeQubit+offset)%16 # this can be used to go clockwise instead
    
    return address

def RunRepetition(bit,d,device,delay,totalRuns):
    
    # set the number of shots to use on the backend
    shots = 8192
    
    # determine whether a simulator is used
    simulator = (device not in ['ibmqx3','ibmqx5'])
    
    # if the simulator is used, we declare the minimum number of qubits required
    if (simulator):
        num = 2*d
    # for the real device there are always 16
    else:
        num = 16
          
    repetitionScript = [] # we create a batch job of totalRuns identical runs
    for run in range(totalRuns):
        
        # set up registers and program
        q = QuantumRegister(num)
        c = ClassicalRegister(num)
        repetitionScript.append( QuantumCircuit(q, c) )

        # now we insert all the quantum gates to be applied
        # a barrier is inserted between each section of the code to prevent the complilation doing things we don't want it to

        # the stored bit is initialized by repeating it across all code qubits same state
        # since qubits are automatically initialized as 0, we just need to do Xs if b=1
        if (bit==1):
            for codeQubit in range(d):
                repetitionScript[run].x( q[GetAddress(codeQubit,0,simulator)] )
            # also do it for the single qubit on the end for comparision
            repetitionScript[run].x( q[GetAddress(d-1,1,simulator)] )

        repetitionScript[run].barrier()

        # if the code is simulated add rotations for error like effects (and a barrier)
        AddDelay(repetitionScript[run],q,num,delay)

        # we then start the syndrome measurements by doing CNOTs between each code qubit and the next ancilla along the line
        for codeQubit in range(d-1):
            AddCnot(repetitionScript[run],q,GetAddress(codeQubit,0,simulator),GetAddress(codeQubit,1,simulator),device)

        repetitionScript[run].barrier()

        # if the code is simulated add rotations for error like effects (and a barrier)
        AddDelay(repetitionScript[run],q,num,delay)

        # next we perform CNOTs between each code qubit and the previous ancilla along the line
        for codeQubit in range(1,d):
            AddCnot(repetitionScript[run],q,GetAddress(codeQubit,0,simulator),GetAddress(codeQubit,-1,simulator),device)

        repetitionScript[run].barrier()

        # if the code is simulated add rotations for error like effects (and a barrier)
        AddDelay(repetitionScript[run],q,num,delay)

        # all qubits are then measured
        for address in range(num):
            repetitionScript[run].measure(q[address], c[address])

        # noise is turned on if simulator is used
        #repetitionScript[run].noise(int(simulator))
            
    # run the job (if the device is available)
    backends = available_backends()
    backend = get_backend(device)
    print('Status of '+device+':',backend.status)

    dataNeeded = True
    while dataNeeded==True:
        if backend.status["available"] is True:
            print("\nThe device is available, so the following job is being submitted.\n")

            print(repetitionScript[1].qasm())
            
            shots_device = 1000
            job = execute(repetitionScript, backend, shots=shots, skip_translation=True)
            
            results = []
            for run in range(totalRuns):
                results.append( job.result().get_counts(repetitionScript[run]) )
            
            dataNeeded = False
        else:
            print("\nThe device is not available, so we will wait for a while.")
            time.sleep(600)
    
    # the raw data states the number of runs for which each outcome occurred
    # we convert this to fractions before output.
    for run in range(totalRuns):
        for key in results[run].keys():
            results[run][key] = results[run][key]/shots
    
    # return the results
    return results

def AddProbToResults(prob,string,results):
    
    if string not in results.keys():
        results[string] = 0
    
    results[string] += prob

def CalculateError(encodedBit,results,table):
    
    # total prob of error will be caculated by looping over all strings
    # some will end up being ignored, so we'll also need to renormalize

    error = 0
    renorm = 1
    
    # all strings from our sample are looped over
    for string in results.keys():

        # we find the probability P(string|encodedBit) from the lookup table
        right = 0
        if string in table[encodedBit].keys():
            right = table[encodedBit][string]
        
        # as is the probability P(string|!encodedBit)
        wrong = 0
        if string in table[(encodedBit+1)%2].keys():
            wrong = table[(encodedBit+1)%2][string]

        # if this is a string for which P(string|!encodedBit)>P(string|encodedBit), the decoding fails
        # the probabilty for this sample is then added to the error
        if (wrong>right):
            error += results[string]
        # if P(string|!encodedBit)=P(string|encodedBit)=0 we have no data to decode, so we should ignore this sample 
        # otherwise if P(string|!encodedBit)=P(string|encodedBit), the decoder randomly chooses between them
        # P(failure|string) is therefore 0.5 in this case
        elif (wrong==right):
            if wrong==0:
                renorm -= results[string]
            else:
                error += 0.5*results[string]
        # otherwise the decoding succeeds, and we don't care about that
    
    if renorm==0:
        error = 1
    else:
        error = error/renorm
            
    return error

def GetData(device,minSize,maxSize,totalRuns,delay):
    
    # loop over code sizes that will fit on the chip (d=3 to d=8)
    for d in range(minSize,maxSize+1):

        print("\n\n**d = " + str(d) + "**")

        # get data for each encoded bit value
        for bit in range(2):

            # run the job and put results in resultsRaw
            results = RunRepetition(bit,d,device,delay,totalRuns)

            delayString = ""
            if delay>0:
                delayString = '_delay='+str(delay)
                
            for run in range(totalRuns):
                f = open('Repetition_Code_Results/'+device+'/results_d=' + str(d) + delayString + '_run=' + str(run) + '_bit=' + str(bit) + '.txt', 'w')
                f.write( str(results[run]) )
                f.close()

def ProcessData(device,encodedBit,minSize,maxSize,totalRuns,delay):
    
    # determine whether a simulator is used
    simulator = (device not in ['ibmqx3','ibmqx5'])
    
    
    # initialize list used to store the calculated means and variances for results from the codes
    codeResults = [[[0]*4 for _ in range(j)] for j in range(minSize,maxSize+1)]
    singleResults = [[[0]*2 for _ in range(16)] for _ in range(minSize,maxSize+1)]
    # singleResults[d-minSize][j][0] is the probability of state 1 for qubit j when used in a code of distance d
    # singleResults[d-minSize][j][1] is the variance for the above
    
    
    # the results will show that the case of partial decoding requires more analysis
    # for this reason we will also output combinedCodeResults, which is all runs of codeResults combined
    # here we initialize list of combined results from the code only case
    combinedResultsCode = [[{} for _ in range(minSize,maxSize+1) ] for _ in range(2)]
    
        
    for d in range(minSize,maxSize+1):
        
        # we loop over code sizes and runs to create the required dictionaries of data:
        # resultsFull, resultsCode and resultsSingle (and well as the things used to make them)
        
        # the results that come fresh from the backend
        resultsVeryRaw = [[{} for _ in range(2)] for run in range(0,totalRuns)]
        resultsRaw = [[{} for _ in range(2)] for run in range(0,totalRuns)]
        # the results from the full code (including ancillas)
        # resultsFull[k] gives results for the effective distance d-k code obtained by ignoring the last k code qubits and ancillas
        resultsFull = [[[{} for _ in range(d)] for _ in range(2)] for run in range(0,totalRuns)]
        # the same but with ancilla results excluded
        resultsCode =  [[[{} for _ in range(d)] for _ in range(2)] for run in range(0,totalRuns)]
        # results each single bit
        resultsSingle = [[[{} for _ in range(16)] for _ in range(2)] for run in range(0,totalRuns)]
        
        
        for run in range(0,totalRuns):
            
            # we get results for both possible encoded bits
            for bit in range(2):
                
                delayString = ""
                if delay>0:
                    delayString = '_delay='+str(delay)
                
                # get results for this run from file
                f = open('Repetition_Code_Results/'+device+'/results_d=' + str(d) + delayString + '_run=' + str(run) + '_bit=' + str(bit) + '.txt')
                resultsVeryRaw[run][bit] = eval(f.read())
                f.close()
                
                
                # loop over all keys in the raw results and look at the ones without strings as values
                # since all such entries should have a bit string as a key, we call it stringVeryRaw
                for stringVeryRaw in resultsVeryRaw[run][bit].keys():
                    if resultsVeryRaw[run][bit][stringVeryRaw] is not str:
                        
                        # create a new dictionary in which each key is padded to a bit string of length 16
                        stringRaw = stringVeryRaw.rjust(16,'0')
                        resultsRaw[run][bit][stringRaw] = resultsVeryRaw[run][bit][stringVeryRaw]


                # now stringRaw only has data in the correct format
                # let's loop over its entries and process stuff
                for stringRaw in resultsRaw[run][bit].keys():

                    # get the prob corresponding to this string
                    probToAdd = resultsRaw[run][bit][stringRaw]

                    # first we deal with resultsFull and resultsCode

                    # loop over all truncated codes relevant for this d
                    for k in range(d):
                        # distance of this truncated code
                        dd = d-k

                        # extract the bit string relevant for resultsFull
                        # from left to right this will alternate between code and ancilla qubits in increasing order
                        stringFull = ''
                        for codeQubit in range(dd): # add bit value for a code qubit...
                            stringFull += stringRaw[15-GetAddress(codeQubit,0,simulator)]
                            if (codeQubit!=(d-1)): #...and then the ancilla next to it (if we haven't reached the end of the code)
                                stringFull += stringRaw[15-GetAddress(codeQubit,1,simulator)]

                        # remove ancilla bits from this to get the string for resultsCode
                        stringCode = ""
                        for n in range(dd):
                            stringCode += stringFull[2*n]
                        AddProbToResults(probToAdd,stringFull,resultsFull[run][bit][k])
                        AddProbToResults(probToAdd,stringCode,resultsCode[run][bit][k])
                    
                    # now we'll do results single

                    # the qubits are listed in the order they are in the code
                    # so for each code qubit
                    for jj in range(8):
                        # loop over it and its neighbour
                        for offset in range(2):
                            stringSingle = stringRaw[15-GetAddress(jj,offset,simulator)]
                            AddProbToResults(probToAdd,stringSingle,resultsSingle[run][bit][2*jj+offset])
                
                # combined this run's resultsCode with the total, using the k=0 values
                for stringCode in resultsCode[run][bit][0].keys():
                    probToAdd = resultsCode[run][bit][0][stringCode]/10
                    AddProbToResults(probToAdd,stringCode,combinedResultsCode[bit][d-minSize])
   
        for run in range(0,totalRuns):

            # initialize list used to store the calculated means and variances for results from the codes
            codeSample = [[0]*2 for _ in range(d)]
            # here
            # codeSample gives the results
            # codeSample[0] gives results for the whole code
            # codeSample[k][0] is the error prob when decoding uses both code and ancilla qubits
            # codeSample[k][1] is the error prob when decoding uses only code qubits
            singleSample = [0]*16
            # singleSample[j] is the probability of state 1 for qubit j when the required bit value is encoded

            # write results in            
            for k in range(d):
                
                # calculate look up tables by averaging over all other runs
                fullTable = [{} for _ in range(2)]
                codeTable = [{} for _ in range(2)]
                for b in range(2):
                    for r in [rr for rr in range(totalRuns) if rr!=run]:
                        for string in resultsFull[r][b][k]:
                            AddProbToResults(resultsFull[r][b][k][string]/(totalRuns-1),string,fullTable[b])
                        for string in resultsCode[r][b][k]:
                            AddProbToResults(resultsCode[r][b][k][string]/(totalRuns-1),string,codeTable[b])
                # then calculate corresponding errors            
                codeSample[k][0] = CalculateError(encodedBit,resultsFull[run][encodedBit][k],fullTable)
                codeSample[k][1] = CalculateError(encodedBit,resultsCode[run][encodedBit][k],codeTable)
                
            for j in range(16):
                if '1' in resultsSingle[run][encodedBit][j].keys():
                    singleSample[j] = resultsSingle[run][encodedBit][j]['1']


            # add results from this run to the overall means and variances
            for k in range(d):
                for l in range(2):
                    codeResults[d-minSize][k][2*l] += codeSample[k][l] / totalRuns # means
                    codeResults[d-minSize][k][2*l+1] += codeSample[k][l]**2 / totalRuns # variances
            for j in range(16):
                singleResults[d-minSize][j][0] += singleSample[j] / totalRuns
                singleResults[d-minSize][j][1] += singleSample[j]**2 / totalRuns

        # finish the variances by subtracting the square of the mean
        for k in range(d):
            for l in range(1,2,4):
                codeResults[d-minSize][k][l] -= codeResults[d-minSize][k][l-1]**2
        for j in range(16):
                singleResults[d-minSize][j][1] -= singleResults[d-minSize][j][0]**2

    
    # return processed results                                                                                                                        
    return codeResults, singleResults, combinedResultsCode

def MakeGraph(X,Y,y,axisLabel,labels=[],legendPos='upper right',verbose=False,log=False,tall=False):
    
    from matplotlib import pyplot as plt
    plt.rcParams.update({'font.size': 30})
    
    markers = ["o","^","h","D","*"]
    
    # if verbose, print the numbers to screen
    if verbose==True:
        print("\nX values")
        print(X)
        for j in range(len(Y)):
            print("\nY values for "+labels[j])
            print(Y[j])
            print("\nError bars")
            print(y[j])
            print("")
    
    # convert the variances of varY into widths of error bars
    for j in range(len(y)):
        for k in range(len(y[j])):
            y[j][k] = math.sqrt(y[j][k]/2)
            
    if tall:
        plt.figure(figsize=(20,20))
    else:
        plt.figure(figsize=(20,10))
    
    
    # add in the series
    for j in range(len(Y)):
        marker = markers[j%len(markers)]
        if labels==[]:
            plt.errorbar(X, Y[j], marker = marker, markersize=20, yerr = y[j], linewidth=5)
        else:
            plt.errorbar(X, Y[j], label=labels[j], marker = marker, markersize=20, yerr = y[j], linewidth=5)
            plt.legend(loc=legendPos)

    
    # label the axes
    plt.xlabel(axisLabel[0])
    plt.ylabel(axisLabel[1])
    
    # make sure X axis is fully labelled
    plt.xticks(X)
    
    
    # logarithms if required
    if log==True:
        plt.yscale('log')

    # make the graph
    plt.show()
    
    plt.rcParams.update(plt.rcParamsDefault)
    

# set device to use
# this also sets the maximum d. We only go up to 6 on the simulator
userInput = input("Do you want results for a real device? (input Y or N) If not, results will be from a simulator. \n").upper()
if (userInput=="Y"):
    device = 'ibmqx3'
else:
    device = 'ibmq_qasm_simulator'

# determine the delay
userInput = input("What value of the delay do you want results for? \n").upper()
delay = 0
try:
    delay = int(userInput)
except:
    pass
    
# determine the code sizes to be considered
userInput = input("What is the minimum size code you wish to consider? (input 3, 4, 5, 6, 7 or 8) \n").upper()
if userInput in ['3','4','5','6','7','8']:
    minSize = int(userInput)
else:
    minSize = 3
    
userInput = input("What is the maximum size code you wish to consider? (input a number no less than the minimum, but no larger than 8) \n").upper()
if userInput in ['3','4','5','6','7','8']:
    maxSize = int(userInput)
else:
    maxSize = 8
    
# determine whether data needs to be taken
if device=='ibmqx5':
    dataAlready = True
else:
    userInput = input("Do you want to process saved data? (Y/N) If not, new data will be obtained. \n").upper()
    if (userInput=="Y"):
        dataAlready = True
    else:
        dataAlready = False

# set number of runs used for stats
# totalRuns is for a repetition code of length d
totalRuns = 10

# if we need data, we get it
if (dataAlready==False):

    # get the required data for the desired number of runs
    GetData(device,minSize,maxSize,totalRuns,delay)

codeResults = [[],[]]
singleResults = [[],[]]
for encodedBit in range(2):
    try:
        codeResults[encodedBit], singleResults[encodedBit], combinedResultsCode = ProcessData(device,encodedBit,minSize,maxSize,totalRuns,delay)
    except Exception as e:
        print(e)

# plot for single qubit data for each code distance
for d in range(minSize,maxSize+1):
    X = range(16)
    Y = []
    y = []
    # a series for each encoded bit
    for encodedBit in range(2):        
        Y.append([singleResults[encodedBit][d-minSize][j][0] for j in range(16)])
        y.append([singleResults[encodedBit][d-minSize][j][1] for j in range(16)])
    # make graph
    print("\n\n***Final state of each qubit for code of distance d = " + str(d) + "***")
    MakeGraph(X,Y,y,['Qubit position in code','Probability of 1'])

for encodedBit in range(2): # separate plots for each encoded bit
    X = range(minSize,maxSize+1)
    Y = []
    y = []
    for dec in range(2): # dec=0 corresponds to full decoding, and 1 to partial
        Y.append([codeResults[encodedBit][d-minSize][0][2*dec+0] for d in range(minSize,maxSize+1)])
        y.append([codeResults[encodedBit][d-minSize][0][2*dec+1] for d in range(minSize,maxSize+1)])
    # minimum error value for the single qubit memory is found and plotted as a comparsion (with max error bars)
    simulator = (device not in ['ibmqx3','ibmqx5'])
    minSingle = min([singleResults[encodedBit][d-minSize][GetAddress(d-1,1,simulator)][0] for d in range(minSize,maxSize+1)])
    maxSingle = max([singleResults[encodedBit][d-minSize][GetAddress(d-1,1,simulator)][1] for d in range(minSize,maxSize+1)])
    Y.append([minSingle]*(maxSize-minSize+1))
    y.append([maxSingle]*(maxSize-minSize+1))
    
    print("\n\n***Encoded " + str(encodedBit) + "***")
    MakeGraph(X,Y,y,['Code distance, d','Error probability, P'],
              labels=['Full decoding','Partial decoding','Single qubit memory'],legendPos='lower left',log=True,verbose=True)
     

for encodedBit in range(2): # separate plots for each encoded bit
    for decoding in ['full','partial']:
        dec = (decoding=='partial') # this is treated as 0 if full and 1 if partial
        X = range(1,maxSize+1)
        Y = []
        y = []
        labels = []
        for d in range(minSize,maxSize+1):# series for each code size
            seriesY = [math.nan]*(maxSize)
            seriesy = [math.nan]*(maxSize)
            for k in range(d):
                seriesY[d-k-1] = codeResults[encodedBit][d-minSize][k][2*dec+0]
                seriesy[d-k-1] = codeResults[encodedBit][d-minSize][k][2*dec+1]
            Y.append(seriesY)
            y.append(seriesy)
            labels.append('d='+str(d))
            
        print("\n\n***Encoded " + str(encodedBit) + " with " + dec*"partial" + (1-dec)*"full" + " decoding***")
        MakeGraph(X,Y,y,['Effective code distance','Logical error probability'],
                  labels=labels,legendPos = 'upper right')   

def MakeModelTables (q,d):
    # outputs an array of two dictionaries for the lookup table for a simple model of a distance d code
    # q[0] is prob of 0->1 noise, and q[1] is prob of 1->0
    # no disinction is made between strings with the same number of errors
    # the prob for all are assigned to a single string, all with 0s on the left and 1s on the right
    
    modelResults = [{},{}]
    bit = ["0","1"]
    for encodedBit in range(2):
        for errors in range(d+1):
            if encodedBit==0:
                string = "0"*(d-errors)+"1"*errors
            else:
                string = "0"*errors+"1"*(d-errors)
            modelResults[encodedBit][string] = scipy.special.binom(d,errors) * q[encodedBit]**errors * (1-q[encodedBit])**(d-errors)
            
    return modelResults

def TotalLogical (p0,d,p1=0): # outputs total logical error prob for a single or two round code
    P0 = CalculateError( encodedBit, MakeModelTables([p0,p0],d)[0], MakeModelTables([p0,p0],d) )
    P1 = CalculateError( encodedBit, MakeModelTables([p0,p0],d)[1], MakeModelTables([p1,p1],d) )
    return P0*(1-P1) + P1*(1-P0)

for encodedBit in range(2): # separate plots for each encoded bit
    
    p = [0]*2 # here is where we'll put p_0 and p_1
    bar = [0]*2
    
    for dec in [1,0]: # dec=0 corresponds to full decoding, and 1 to partial
        
        # get the results we want to fit to
        realResults = [codeResults[encodedBit][d-minSize][0][2*dec] for d in range(minSize,maxSize+1)]
         
        # search possible values intul we find a minimum (assumed to be global)
        # first we do the partial decoding to get p (which is stored in p[0])
        # then the full decoding to get p[1]=p_1
        minimizing = True
        delta = 0.001
        q = delta
        diff =[math.inf,math.inf]
        while minimizing:
            
            q += delta # set new q
            diff[0] = diff[1] # copy diff value for last p
            
            # calculate diff for new q
            diff[1] = 0
            for d in range(minSize,maxSize+1):
                if dec==1:
                    Q = TotalLogical(q,d)
                else:
                    Q = TotalLogical(p[0]-q,d,p1=q)
                diff[1] += ( math.log( realResults[d-minSize] ) - math.log( Q ) )**2
            # see if a minimum has been found
            minimizing = ( diff[0]>diff[1] )
            
        # go back a step on p to get pSum
        p[1-dec] = q - delta
        
        # get diff per qubit
        bar[1-dec] = math.exp(math.sqrt( diff[0]/(maxSize-minSize+1) ))
    
    p[0] = p[0] - p[1] # put p_0 in p[0] (instead of p)
        
    print("\n\n***Encoded " + str(encodedBit) + "***\n" )
        
    for j in [0,1]:
        print("   p_"+str(j)+" = " + str(p[j]) + " with fit values typically differing by a factor of " + str(bar[j]) + "\n")

        
    plottedMinSize = max(4,minSize) # we won't print results for d=3 for clarity
    X = range(plottedMinSize,maxSize+1)
    Y = []
    y = []
    # original plots
    for dec in range(2): # dec=0 corresponds to full decoding, and 1 to partial
        # results from the device
        Y.append([codeResults[encodedBit][d-minSize][0][2*dec+0] for d in range(plottedMinSize,maxSize+1)])
        y.append([codeResults[encodedBit][d-minSize][0][2*dec+1] for d in range(plottedMinSize,maxSize+1)])
    # fit lines
    for dec in range(2):
        if dec==1:
            Y.append([TotalLogical(p[0]+p[1],d) for d in range(plottedMinSize,maxSize+1)])
        else:
            Y.append([TotalLogical(p[0],d,p1=p[1]) for d in range(plottedMinSize,maxSize+1)])
        y.append([0]*(maxSize-plottedMinSize+1))

    MakeGraph(X,Y,y,['Code distance, d','Error probability, P'],
              labels=['Full decoding','Partial decoding','Full decoding (model)','Partial decoding (model)'],legendPos='lower left',log=True)
     

# for each code distance and each encoded bit value, we'll create a list of the probabilities for each possible number of errors
# list is initialized with zeros
errorNum = [[[0]*(d+1) for d in range(minSize,maxSize+1)] for _ in range(2)]

for d in range(minSize,maxSize+1):
        for bit in range(2):
            # for each code distance and each encoded bit value we look at all possible result strings
            for string in combinedResultsCode[bit][d-minSize]:
                # count the number of errors in each string
                num = 0
                for j in range(d):
                    num += ( int( string[j] , 2 ) + bit )%2
                # add prob to corresponding number of errors
                errorNum[bit][d-minSize][num] += combinedResultsCode[bit][d-minSize][string]
        

        # the we make a graph for each, and print a title
        Y0 = [y if y>0 else math.nan for y in errorNum[0][d-minSize]]
        Y1 = [y if y>0 else math.nan for y in errorNum[1][d-minSize]]
        print("\n\n***Probability of errors on code qubits for d = " + str(d) + "***")
        MakeGraph(range(d+1),[Y0,Y1],[[0]*(d+1)]*2,['Number of code qubit errors','Probability (log base 10)'],
                  labels=['Encoded 0','Encoded 1'],legendPos='upper right',log=True)

        # actually, we make two graphs. This one plots the number of 1s rather than errors, and so the plot for encoded 1 is inverted
        # Y0 in this graph is as before
        Y1 = Y1[::-1] # but Y1 has its order inverted
        print("\n\n***Probability for number of 1s in code qubit result for d = " + str(d) + "***")
        MakeGraph(range(d+1),[Y0,Y1],[[0]*(d+1)]*2,['Number of 1s in code qubit result','Probability (log base 10)'],
                  labels=['Encoded 0','Encoded 1'],legendPos='center right',log=True)    