How to extend the number of Symfit parameters arbitrarily in Python

Question

I have the following optimization code that is a parameterized by a variable n.

from symfit import parameters, Eq, Ge, Fit
import numpy as np
n = 3
xdata = np.sort(np.random.choice(range(1, 4*n), n)) # Make fake data
print(xdata)
p1, p2, p3 = parameters('p1, p2, p3')
model = p1*p2*p3
constraints = [
    Eq(xdata[0]*p1+(xdata[1]-xdata[0])*p2+(xdata[2]-xdata[1])*p3, 1),
    Ge(p1, p2),
    Ge(p2, p3),
    Ge(p3, 0)
    ]

fit = Fit(- model, constraints=constraints)
fit_result = fit.execute()
print(fit_result)

I would like to use it for much larger values of n but I don't know how to change the lines

 p1, p2, p3 = parameters('p1, p2, p3')
 model = p1*p2*p3

and the constraints to cope with an arbitrarily large n.

The code is using the symfit library. That link shows an example of how parameters is used and a link to the documentation.

How can one do this?

Answer 1

Numpy interacts really well with the symfit library. All of the operations you are trying to generalize are fairly trivial when using it.

---

Setup

n = 3
_data = np.sort(np.random.choice(np.arange(1, 4 * n), n))

---

1. String formatting

   You can dynamically create a tuple of parameters using a simple iterator and str.join, which you can then pass into the parameters constructor to get a tuple of your parameters.

params = parameters(', '.join(f'p{i}' for i in range(1, n+1)))
                ^^
# p1, p2, p3 = parameters('p1, p2, p3')

2. np.prod

   This operation is quite straight forward. np.prod computes the:

   product of array elements over a given axis

   Which when applies to a tuple of symfit parameters, produces your desired p1*p2*...pn

model = np.prod(params)
        ^^
# model = p1*p2*p3

3. np.concatenate + np.diff

   Probably the most complicated line to generalize, but still not too hard to understand. You want to multiply the differences of consecutive elements in the data array by your parameters, and sum the results. Since the first element will not have a difference with a previous element, you can use np.concatenate to add it back in.

u = np.concatenate((_data[:1], np.diff(_data)))
c1 = Eq(np.sum(params * u), 1)
              ^^
# Eq(xdata[0]*p1+(xdata[1]-xdata[0])*p2+(xdata[2]-xdata[1])*p3, 1)

4. np.column_stack

   You want a rolling view of your parameters as constraints: p1-p2, p2-p3, ...pn, 0. This just stacks a one-off tuple with a zero padded on with the original tuple of parameters, then uses a list comprehension to unpack into your Ge constructors.

ges = [Ge(*group) for group in np.column_stack((params, params[1:] + (0,)))]

5. Fit!

   I didn't change anything here!

constraints = [c1, *ges]
fit = Fit(- model, constraints=constraints)
fit_result = fit.execute()