opulation size: 1000
Calculated number of qubits: 10
Valid population fitnesses: [0.74608824 0.85743972 ........... 0.73538867 0.50542962]
Running Grover search with 10 qubits
Initial max fitness threshold: 3320

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 123
    121 # Run Grover’s search to find the optimal individual
    122 print(f"Valid population fitnesses: {valid_fitnesses}")  # Debugging print statement
--> 123 best_index, optimal_solution = run_grover_search(num_qubits, valid_fitnesses, num_iterations)
    125 # Display the optimal individual and its index
    126 print(f"Optimal solution (best input parameters for individual {best_index}): {optimal_solution}")

Cell In[3], line 100, in run_grover_search(num_qubits, valid_fitnesses, num_iterations)
     97 # Apply fitness evaluation unitary and Grover iterations
     98 for _ in range(num_iterations):
     99     # Apply the fitness evaluation step with modified oracle
--> 100     oracle = create_two_complement_oracle(valid_fitnesses, num_qubits, max_fitness)
    101     if oracle.size() > 0:  # Ensure oracle is non-trivial
    102         grover_op = GroverOperator(oracle)

Cell In[3], line 60, in create_two_complement_oracle(valid_fitnesses, total_qubits, current_max)
     57 threshold_bin = format(current_max + 1, f'0{total_qubits}b')
     59 if len(threshold_bin) > total_qubits:
---> 60     raise ValueError(f"Threshold exceeds the available number of qubits: {len(threshold_bin)} vs {total_qubits}")
     62 # Subtract (l+1) using two's complement arithmetic
     63 for i, bit in enumerate(threshold_bin):

ValueError: Threshold exceeds the available number of qubits: 12 vs 10


My code:
import numpy as np
from qiskit import QuantumCircuit, ClassicalRegister, transpile
from qiskit.quantum_info import Statevector
from qiskit.circuit.library import GroverOperator
from math import ceil, log2

# Parameters for population
num_params = 12  # Number of scalar variables per individual
population_size = 1000  # Number of individuals in the population
num_iterations = 3  # Grover's iterations

# Random population of 1000 individuals, each having 12 scalar variables
population = np.random.rand(population_size, num_params)
print(f"Population size: {len(population)}")  # Debugging print statement

# Random fitness values for testing (some valid, some invalid, denoted by NaN)
population_fitnesses = np.random.rand(population_size)
invalid_indices = np.random.choice(population_size, size=int(0.2 * population_size), replace=False)  # 20% invalid
population_fitnesses[invalid_indices] = np.nan  # Marking invalid fitnesses as NaN

# Filter valid fitness values
valid_indices = np.where(~np.isnan(population_fitnesses))[0]
valid_fitnesses = population_fitnesses[valid_indices]

# Number of qubits required (based on the valid population size)
num_qubits = ceil(log2(len(valid_fitnesses)))
print(f"Calculated number of qubits: {num_qubits}")  # Debugging print statement

# Initialize quantum population with superposition
def initialize_population(num_qubits):
    qc = QuantumCircuit(num_qubits)
    qc.h(range(num_qubits))  # Hadamard gate to create equal superposition of all states
    return qc

# Apply the fitness function to update the fitness register
def fitness_evaluation(qc, valid_fitnesses, total_qubits):
    oracle_circuit = QuantumCircuit(total_qubits + 1)  # Fitness register gets extra qubit for validity
    
    for i, fitness in enumerate(valid_fitnesses):
        bin_fitness = format(int(fitness), f'0{total_qubits}b')  # Fitness encoded as binary
        
        # Apply X-gates to flip bits where necessary
        for j, bit in enumerate(bin_fitness):
            if bit == '0':
                oracle_circuit.x(j)

        # Mark invalid fitnesses (i.e., lower bound for valid range)
        oracle_circuit.x(total_qubits)  # Last bit as '0' for invalid

    return oracle_circuit

# Oracle for marking individuals with fitness greater than current_max using two's complement addition
def create_two_complement_oracle(valid_fitnesses, total_qubits, current_max):
    oracle_circuit = QuantumCircuit(total_qubits)
    
    # Ensure that the current_max threshold is represented using the correct number of bits
    threshold_bin = format(current_max + 1, f'0{total_qubits}b')
    
    if len(threshold_bin) > total_qubits:
        raise ValueError(f"Threshold exceeds the available number of qubits: {len(threshold_bin)} vs {total_qubits}")

    # Subtract (l+1) using two's complement arithmetic
    for i, bit in enumerate(threshold_bin):
        if bit == '1' and i < total_qubits:  # Ensure the index doesn't exceed the qubit count
            oracle_circuit.x(i)  # Apply X gate to implement subtraction in two's complement

    # Oracle marks the states with fitness greater than the threshold
    for i, bit in enumerate(threshold_bin):
        if bit == '0' and i < total_qubits:
            oracle_circuit.x(i)  # Flip the bits of states below the threshold

    # Restore the original fitness values by adding (l+1) back
    for i, bit in enumerate(threshold_bin):
        if bit == '1' and i < total_qubits:
            oracle_circuit.x(i)  # Undo the subtraction

    return oracle_circuit

# Grover's search for optimizing parameters
def run_grover_search(num_qubits, valid_fitnesses, num_iterations):
    required_qubits = ceil(log2(len(valid_fitnesses)))
    if num_qubits < required_qubits:
        raise ValueError(f"Not enough qubits to represent valid population indices. Required: {required_qubits}, but got {num_qubits}.")
    
    print(f"Running Grover search with {num_qubits} qubits")  # Debugging print statement

    # Initialize the quantum circuit
    qc = QuantumCircuit(num_qubits)

    # Create a superposition of all valid individuals
    qc.compose(initialize_population(num_qubits), inplace=True)

    # Initialize the maximum fitness threshold randomly between 2^(M+1) and 2^(M+2)-1
    max_fitness = np.random.randint(2**(num_qubits+1), 2**(num_qubits+2))
    print(f"Initial max fitness threshold: {max_fitness}")

    # Apply fitness evaluation unitary and Grover iterations
    for _ in range(num_iterations):
        # Apply the fitness evaluation step with modified oracle
        oracle = create_two_complement_oracle(valid_fitnesses, num_qubits, max_fitness)
        if oracle.size() > 0:  # Ensure oracle is non-trivial
            grover_op = GroverOperator(oracle)
            qc.compose(grover_op, inplace=True)
        
        # Apply Grover's algorithm iterations for amplification
        max_fitness_index = fitness_evaluation(qc, valid_fitnesses, num_qubits)
        
        # Update max_fitness if a better solution is found
        if valid_fitnesses[max_fitness_index] > max_fitness:
            max_fitness = valid_fitnesses[max_fitness_index]
            print(f"Updated max fitness threshold: {max_fitness}")

    # Measure the result
    statevector = Statevector.from_instruction(qc)
    probabilities = statevector.probabilities_dict()
    best_state = max(probabilities, key=probabilities.get)
    best_index = int(best_state, 2) % len(valid_fitnesses)

    return valid_indices[best_index], population[valid_indices[best_index]]

# Run Grover’s search to find the optimal individual
print(f"Valid population fitnesses: {valid_fitnesses}")  # Debugging print statement
best_index, optimal_solution = run_grover_search(num_qubits, valid_fitnesses, num_iterations)

# Display the optimal individual and its index
print(f"Optimal solution (best input parameters for individual {best_index}): {optimal_solution}")