I am using nlopt in Python, I’m taking some values in a matrix, defining a function in as many variables as the sum of its two dimensions, setting up some constraints, and optimizing. This helps me find the parameters that maximize a graph entropy expression.

Adding a variable number of constraints based on this matrix is proving difficult, I keep getting “index 0 is out of bounds” errors:

    result[count] = vars[i] * vars[j] - 1  # vars[i] * vars[j] < 1
    ~~~~~~^^^^^^^
    IndexError: index 0 is out of bounds for axis 0 with size 0
   

when I try and add the constraints. Essentially, what is wrong with this method of adding constraints?

import numpy as np
import sympy as sp
import ast
import nlopt

Qgreen = np.array([[1, 2, 3], 
                  [4, 5, 6], 
                  [7, 8, 9]])

len1, len2 = Qgreen.shape

# Calculate strength along the rows
strength1 = [np.sum(Qgreen[i, :]) for i in range(len1)]

# Calculate strength along the columns
strength2 = [np.sum(Qgreen[:, i]) for i in range(len2)]

# Combine the strengths into one list
str_combined = strength1 + strength2

# Define symbolic variables
vars = sp.symbols(f'x1:{len(str_combined)+1}')  # x1, x2, ..., x6

# Create the function using symbolic expressions
f_sum = sum(s * sp.log(v) for s, v in zip(str_combined, vars))

# Flatten the product terms and their logarithms
product_terms = []
for i in range(len1):
    for j in range(len1, len1 + len2):
        product_terms.append(vars[i] * vars[j])

# Adding the log(1 - product_terms) part of the expression
f_product = sum(sp.log(1 - term) for term in product_terms)

# Final function combining both parts
f = f_sum + f_product

# Create a lambda function for fast evaluation
f_lambdified = sp.lambdify(vars, f, "numpy")

# Define the function for nlopt (minimization)
def nlopt_function(vars, grad):
    if grad.size > 0:
        # Compute the gradient if needed
        grad[:] = np.array([sp.diff(f, var).evalf(subs={v: vars[i] for i, v in enumerate(vars)}) for i, var in enumerate(vars)])
    return f_lambdified(*vars)

def constraints(vars, result):
    count = 0
    num_vars = len(vars)  # Total number of variables
    print(num_vars)
    for i in range(num_vars):
        for j in range(num_vars):
                if i != j:  # Avoid self-products
                    result[count] = vars[i] * vars[j] - 1  # The constraints are vars[i] * vars[j] < 1
                    count += 1

opt = nlopt.opt(nlopt.LN_COBYLA, len(vars))  # Use COBYLA which supports constraints
opt.set_max_objective(nlopt_function)
opt.add_inequality_constraint(constraints, 1e-8)
initial_guess = np.random.uniform(0, 1, len(vars))
optimized_vars = opt.optimize(initial_guess)
max_val = opt.last_optimum_value()
print("Optimized variables:", optimized_vars)
print("Maximum value:", max_val)

I’m trying various combinations of constraint adding, but keep getting the same index out of bounds error.