Entropy Computing via QCI (Python)

# eqc-models 0.20.2 omits its optional direct-hardware client on Python 3.11+.
# Suppress only that exact import-time log; cloud clients remain available.
import logging


class OptionalDirectSolverNotice(logging.Filter):
    def filter(self, record):
        return not (
            record.name == "root"
            and record.levelno == logging.WARNING
            and record.getMessage() == "eqc-direct package not available"
            and record.pathname.replace("\\", "/").endswith("/eqc_models/solvers/eqcdirect.py")
        )


_direct_notice = OptionalDirectSolverNotice()
logging.getLogger().addFilter(_direct_notice)
try:
    import eqc_models
    from eqc_models.base import ConstrainedPolynomialModel, PolynomialModel, QuadraticModel
    from eqc_models.base.operators import Polynomial
    from eqc_models.solvers import Dirac3IntegerCloudSolver, Dirac3ContinuousCloudSolver
finally:
    logging.getLogger().removeFilter(_direct_notice)

# Import Matplotlib to generate plots
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker

# Import numpy and scipy for certain numerical calculations below
import numpy as np
import pyomo.environ as pyo
import itertools
import io
import os
import warnings
from contextlib import redirect_stderr, redirect_stdout

warnings.filterwarnings(
    "ignore",
    message="Cannot evaluate objective value in results.*",
    category=UserWarning,
)


def solve_without_provider_identifiers(solver, model, **kwargs):
    """Submit a QCI job without publishing provider job or file identifiers."""
    with redirect_stdout(io.StringIO()), redirect_stderr(io.StringIO()):
        response = solver.solve(model, **kwargs)
    print("QCI job completed; provider identifiers are omitted from notebook output.")
    return response


class ScalarConstrainedPolynomialModel(ConstrainedPolynomialModel):
    @staticmethod
    def _scalar(value):
        return float(np.asarray(value).reshape(-1)[0])

    def evaluateObjective(self, solution):
        return self._scalar(super().evaluateObjective(solution))

    def evaluatePenalties(self, solution, include_offset=False):
        coefficients, indices = self.penalties
        value = self.offset if include_offset else 0
        value += self._scalar(
            Polynomial(coefficients, indices).evaluate(
                np.array(solution, dtype=np.float64)
            )
        )
        return float(value)

api_url = "https://api.qci-prod.com"
enable_qci_cloud = os.environ.get("QUBONOTEBOOKS_QCI_ENABLE_CLOUD", "0") == "1"
api_token = ""
if enable_qci_cloud:
    api_token = os.environ.get("QCI_TOKEN", "")
    if IN_COLAB and not api_token:
        from google.colab import userdata

        api_token = userdata.get("QCI_TOKEN") or ""
    if not api_token:
        raise RuntimeError("Set QCI_TOKEN before opting into QCi cloud execution.")
else:
    print("QCi cloud disabled; all model construction and local checks will run.")

# 1. Create a concrete Pyomo model
m = pyo.ConcreteModel(name="Simple_Quadratic_Program")

# 2. Define the decision variables with their bounds
# The variables x_1, x_2, and x_3 are defined with bounds [0, 10].
m.x = pyo.Var([1, 2, 3], domain=pyo.NonNegativeReals, bounds=(0, 10))

# 3. Define the objective function
# The objective is to minimize the quadratic expression.
m.obj = pyo.Objective(
    expr=3*m.x[1]**2 + 2*m.x[2]**2 + m.x[3]**2, 
    sense=pyo.minimize
)

# 4. Define the linear constraint
# The sum of the variables must equal 10.
m.c1 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 10)

# 5. Create a solver instance and solve the model
# HiGHS solves this convex quadratic program without an external executable.
solver = pyo.SolverFactory('highs')
if not solver.available(exception_flag=False):
    raise RuntimeError("HiGHS not found. Install the locked QCi environment above.")
results = solver.solve(m, tee=False)
pyo.assert_optimal_termination(results)

# All variables and objective values in these examples are dimensionless.
# The equality-constrained analytic solution satisfies 6*x1 = 4*x2 = 2*x3.
continuous_reference = np.array([pyo.value(m.x[i]) for i in m.x])
np.testing.assert_allclose(continuous_reference, [20/11, 30/11, 60/11], rtol=1e-6)
np.testing.assert_allclose(pyo.value(m.obj), 600/11, rtol=1e-6)

# 6. Display the optimization results
print("\n" + "="*30)
print("-- 📊 Optimization Results --")
print(f"Solver Status: {results.solver.status}")
print(f"Termination Condition: {results.solver.termination_condition}")
print(f"Objective Value (E): {pyo.value(m.obj):.4f}")
print("\n-- 🏗️ Variable Values --")
for i in m.x:
    print(f"x[{i}] = {pyo.value(m.x[i]):.4f}")
print("="*30)

R = 10
coefficients = [3, 2, 1]
indices = [(1,1), (2,2), (3,3)] # Create a polynomial model
model = PolynomialModel(coefficients, indices)
model.upper_bound = R*np.ones((3,)) # Bounds for the variables
# Verify the provider's local polynomial representation against the reference.
np.testing.assert_allclose(np.asarray(model.evaluate(continuous_reference)).item(), 600/11, rtol=1e-6)

continuous_response = None
if enable_qci_cloud:
    solver = Dirac3ContinuousCloudSolver(url=api_url, api_token=api_token)
    continuous_response = solve_without_provider_identifiers(
        solver,
        model,
        sum_constraint=10,
        relaxation_schedule=1,
        num_samples=20,
    ) # sum constraint is set to 10, 20 samples are collected

    # Print the results
    print(continuous_response["results"]["solutions"])
    print(continuous_response["results"]["energies"])
    print(continuous_response["results"]["counts"])
else:
    print("Skipped QCi cloud submission.")

if continuous_response is not None:
    print({
        key: continuous_response["results"][key]
        for key in ("solutions", "energies", "counts")
    })

if continuous_response is not None:
    # Extract the results from the continuous_response
    results = continuous_response["results"]
    energies = results["energies"]
    solutions = results["solutions"]

    # Find the index of the best solution (minimum energy)
    best_index = np.argmin(energies)

    # Get the best solution and its energy
    best_energy = energies[best_index]
    best_solution = solutions[best_index]

    # Print the results in a clear format
    print("-- 📊 Best Solution Found --")
    print(f"Optimal Objective (Energy): {best_energy:.4f}")
    print(f"Variable Values: {np.round(best_solution, 4)}")

# Enumerate the four binary assignments directly.
# For a 2-variable QUBO this is clearer and avoids an external MINLP solver.
def simple_qubo_energy(x1, x2):
    return -0.75 * x1 - 0.25 * x2 + 2 * x1 * x2 + 0.5 * x1

assignments = []
for x1, x2 in itertools.product([0, 1], repeat=2):
    assignments.append({"x": (x1, x2), "energy": simple_qubo_energy(x1, x2)})

best_assignment = min(assignments, key=lambda item: item["energy"])

print("\n" + "="*30)
print("-- Enumeration Results --")
for assignment in assignments:
    x1, x2 = assignment["x"]
    print(f"x[1] = {x1}, x[2] = {x2}, objective = {assignment['energy']:.4f}")
print(f"Optimal Objective: {best_assignment['energy']:.4f}")
print(f"Variable Values: x = {list(best_assignment['x'])}")
print("="*30)

coefficients = [-0.75, -0.25, 2, 0.5]  # Coefficients for the polynomial model
indices = [(1, 1), (2, 2), (1, 2), (0, 1)]  # Create a polynomial model
model = PolynomialModel(coefficients, indices)

# add upper bounds for the variables
model.upper_bound = np.ones((2,))  # Bounds for the variables
# Compare all four binary energies, including both degenerate minima.
binary_candidates = np.array(list(itertools.product([0, 1], repeat=2)))
binary_energies = np.asarray([model.evaluate(x) for x in binary_candidates]).reshape(-1)
np.testing.assert_allclose(binary_energies, [0, -0.25, -0.25, 1.5])

binary_response = None
if enable_qci_cloud:
    solver = Dirac3IntegerCloudSolver(url=api_url, api_token=api_token)
    binary_response = solve_without_provider_identifiers(
        solver,
        model,
        relaxation_schedule=1,
        num_samples=100,
    )

    # Print the results
    print(binary_response["results"]["solutions"])
    print(binary_response["results"]["energies"])
    print(binary_response["results"]["counts"])
else:
    print("Skipped QCi cloud submission.")

def plot_solution_distribution(response):
    """
    Processes a solver response and plots the distribution of unique solutions.

    Args:
        response (dict): The response dictionary from the Dirac-3 solver,
                         containing 'solutions' and 'counts' in its 'results' key.
    """
    # --- Data Processing ---
    # Extract the raw results from the response dictionary
    results = response.get('results', {})
    raw_solutions = results.get('solutions', [])
    raw_counts = results.get('counts', [])

    if not raw_solutions or not raw_counts:
        print("Response object does not contain valid solutions or counts.")
        return

    # Aggregate counts for each unique solution
    solution_counts = {}
    for solution, count in zip(raw_solutions, raw_counts):
        # Convert list to tuple to use it as a dictionary key
        solution_tuple = tuple(solution)
        solution_counts[solution_tuple] = solution_counts.get(solution_tuple, 0) + count

    # Prepare data for plotting
    # Create string labels for the x-axis from the solution tuples
    plot_solutions_str = [str(list(s)) for s in solution_counts.keys()]
    plot_counts = list(solution_counts.values())

    # --- Plotting ---
    plt.style.use('seaborn-v0_8-whitegrid')
    fig, ax = plt.subplots(figsize=(8, 7)) # Increased height for rotated labels

    # Create the bar plot
    bars = ax.bar(plot_solutions_str, plot_counts, color='#4A90E2', width=0.5)

    # Add labels and title for clarity
    ax.set_xlabel("Solution", fontsize=12, labelpad=15)
    ax.set_ylabel("Counts (Frequency)", fontsize=12)
    ax.set_title("Distribution of Solutions", fontsize=14, fontweight='bold')
    
    # Rotate x-axis labels to prevent overlap
    plt.xticks(rotation=90)
    
    # Ensure the y-axis only displays integer values
    ax.yaxis.set_major_locator(mticker.MaxNLocator(integer=True))

    # Add count labels on top of each bar
    for bar in bars:
        yval = bar.get_height()
        ax.text(bar.get_x() + bar.get_width()/2.0, yval + 0.1, int(yval), ha='center', va='bottom', fontsize=11)

    # Set y-axis to start at 0 and give some space at the top
    if plot_counts:
        ax.set_ylim(0, max(plot_counts) * 1.15)

    # Adjust layout to make room for rotated labels
    plt.tight_layout()
    plt.show()

if binary_response is not None:

    plot_solution_distribution(binary_response)

A = np.array([[1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1],
            [0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1],
            [0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1]])
b = np.array([1, 1, 1])
c = np.array([2, 4, 4, 4, 4, 4, 5, 4, 5, 6, 5])
num_variables = A.shape[1]

# Build a Pyomo model for the constrained linear integer program
m = pyo.ConcreteModel(name="Constrained_Linear_Integer_Program")

# Define the set of variable indices
m.I = pyo.RangeSet(0, num_variables - 1)
# Define the set of constraint indices
m.J = pyo.RangeSet(0, A.shape[0] - 1)

# Define the 11 binary variables
m.x = pyo.Var(m.I, domain=pyo.Binary)

# Define the Objective Function
# The objective is the original linear expression.
m.obj = pyo.Objective(
    expr=sum(c[i] * m.x[i] for i in m.I),
    sense=pyo.minimize
)

# Define the constraints
# We define a set of constraints Ax = b directly.
def ax_constraint_rule(model, j):
    # For each row j, sum(A[j,i] * x[i]) must equal b[j]
    return sum(A[j, i] * model.x[i] for i in model.I) == b[j]

m.constraints = pyo.Constraint(m.J, rule=ax_constraint_rule)

# Solve the model
# We use a Mixed-Integer Linear Programming (MILP) solver.
solver = pyo.SolverFactory('highs')
if not solver.available(exception_flag=False):
    raise RuntimeError("HiGHS not found. Install the locked QCi environment above.")
results = solver.solve(m, tee=False)
pyo.assert_optimal_termination(results)

# Exhaustive enumeration is inexpensive for 11 dimensionless binary variables.
local_candidates = np.array(list(itertools.product([0, 1], repeat=num_variables)))
feasible_mask = np.all(local_candidates @ A.T == b, axis=1)
feasible_costs = local_candidates[feasible_mask] @ c
local_optimum = int(feasible_costs.min())
assert local_optimum == 5
np.testing.assert_allclose(pyo.value(m.obj), local_optimum)
np.testing.assert_allclose(A @ np.array([pyo.value(m.x[i]) for i in m.I]), b)

# Display the optimization results
print("\n" + "="*40)
print("-- 📊 Optimization Results --")
print(f"Solver Status: {results.solver.status}")
print(f"Termination Condition: {results.solver.termination_condition}")
print(f"Objective Value: {pyo.value(m.obj):.4f}")

print("\n-- 🏗️ Variable Values --")
solution_vector = [int(pyo.value(m.x[i])) for i in m.I]
print(f"x = {solution_vector}")

epsilon = 1 # choose rho strictly above the sufficient penalty bound
rho = np.sum(np.abs(c)) + epsilon # penalty multiplier

# Construct the quadratic Hamiltonian J
J = rho*np.matmul(A.T,A)

# Construct the linear Hamiltonian C
C = c.copy()  
C -= rho*2*(np.matmul(b.T,A))

# Construct the offset term cQ
cQ = rho*np.matmul(b.T,b)

# --- Formatted Printing ---
print("\n" + "="*60)
print("--- QUBO Formulation Parameters ---")
# Set numpy print options for better readability
np.set_printoptions(precision=2, suppress=True)

print("\nQuadratic Hamiltonian (J):\n", J)
print("\n" + "-"*60)
print("\nLinear Hamiltonian (C):\n", C)
print("\n" + "-"*60)
print(f"\nOffset Term (cQ): {cQ:.2f}")
print("\n" + "="*60)

model = QuadraticModel(C, J)
model.upper_bound = np.ones((11, )) # Set the upper bound for the variables

# Match every energy to the original objective plus squared residual penalties.
manual_energies = np.array([model.evaluate(x) + cQ for x in local_candidates])
reference_energies = local_candidates @ c + rho * np.sum((local_candidates @ A.T - b)**2, axis=1)
np.testing.assert_allclose(manual_energies, reference_energies)
local_minimizers = local_candidates[manual_energies == manual_energies.min()]
assert manual_energies.min() == local_optimum
assert np.all(local_minimizers @ A.T == b)
print(f"Validated {len(local_candidates)} QUBO energies; local optimum = {local_optimum}.")

qubo_response = None
if enable_qci_cloud:
    solver = Dirac3IntegerCloudSolver(url=api_url, api_token=api_token)
    Q = model.qubo.Q

    qubo_response = solve_without_provider_identifiers(solver, model, num_samples=10) # qubomodel

    solution = model.decode(np.array(qubo_response["results"]["solutions"][0]), "qubo")

    true_solution = model.evaluate(solution) + cQ # add the offset term to get the true objective value

    print(f"Decoded Solution (x): {solution}")
    print(f"Objective Value (E): {true_solution}")
else:
    print("Skipped QCi cloud submission; QUBO checked locally.")

if qubo_response is not None:
    # Plot the distribution of unique solutions
    plot_solution_distribution(qubo_response)

coefficients = c
num_vars = 11
indices = [[0, i + 1] for i in range(num_vars)]

lhs = A
rhs = b

# ConstrainedPolynomialModel wraps the QUBO with explicit constraints.
constraint_model = ScalarConstrainedPolynomialModel(coefficients, indices, lhs, rhs)
constraint_model.upper_bound = np.ones(num_vars, ) # Set the upper bound for the variables
# Set the penalty multiplier for the constraints
# This is a crucial step to ensure the constraints are satisfied in the optimization process.
constraint_model.penalty_multiplier = np.sum(np.abs(c)) + 1  # Set the penalty multiplier
print(constraint_model.offset) # Offset term generated by the model

# Check the independent constrained representation, including its offset.
constrained_energies = np.asarray([
    constraint_model.evaluate(x) for x in local_candidates
]).reshape(-1) + constraint_model.offset * constraint_model.penalty_multiplier
np.testing.assert_allclose(constrained_energies, reference_energies)
for candidate in local_candidates:
    np.testing.assert_allclose(constraint_model.evaluateObjective(candidate), c @ candidate)
    np.testing.assert_allclose(
        constraint_model.evaluatePenalties(candidate, include_offset=True),
        np.sum((A @ candidate - b)**2),
    )
print("Validated constrained objective, penalties, and offset for all binary assignments.")

constrained_response = None
if enable_qci_cloud:
    solver = Dirac3IntegerCloudSolver(url=api_url, api_token=api_token)
    constrained_response = solve_without_provider_identifiers(
        solver,
        constraint_model,
        num_samples=10,
    )

else:
    print("Skipped QCi cloud submission; constrained model checked locally.")

if constrained_response is not None:
    # plot the distribution of unique solutions
    plot_solution_distribution(constrained_response)

    # Extract the best solution and its objective value
    best_solution = constrained_response["results"]["solutions"][0]
    best_objective_value = constrained_response["results"]["energies"][0] + constraint_model.offset * constraint_model.penalty_multiplier 
    # Add the offset term to get the true objective value

    print(f"Best Solution: {best_solution}")
    print(f"Best Objective Value: {best_objective_value}")

# =============================================================================
# EXERCISE 1: Adjust the sum constraint
# =============================================================================
# Change the sum constraint in the small QCi example and predict how the feasible solution set changes before submitting or enumerating the model.
#
# Hint: For credential-free validation, enumerate a tiny version locally.
#
# Your code here:
# =============================================================================

# =============================================================================
# EXERCISE 2: Validate a small QUBO locally
# =============================================================================
# Before submitting a two-variable model, enumerate all binary assignments and compute the energy by hand or with a short helper.
#
# Hint: There are only four assignments.
#
# Your code here:
# =============================================================================

# =============================================================================
# EXERCISE 3: Plan for missing credentials
# =============================================================================
# Write down which cells can run without QCi credentials and which must be skipped. Add a fallback analysis that still checks model construction.
#
# Hint: Separate model-building validation from cloud submission.
#
# Your code here:
# =============================================================================