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()

Answer 2

The parameters string needs to be computed in a dynamic manner from n

paramstr = ', '.join(['p{}'.format(i) for i in range(1, n)])
# when n=1, paramstr = 'p1, p2, p3'

Use paramstr as the argument to parameters function.

paramvals = parameters(paramstr)

model can be refactored by reducing paramvals over its product.

from functools import reduce

model = reduce(lambda x, y: x * y, paramvals, 1)

Now to the interesting bit! The constraints can be refactored as:

eqs = xdata[0] * paramvals[0] + sum(
    (xdata[i] - xdata[i-1]) * paramvals[i]
    for i in range(1, n)
)

ges = [
    Ge(paramvals[i-1], paramvals[i])
    for i in range(1, n)
]

ges.append(
    Ge(paramvals[-1], 0)
)

constraints = [
    Eq(eqs, 1),
    *ges
]

Answer 3

I don't know anything about Symfit, but if you're simply trying to generalize the above code to arbitrary N, then:

• You can arbitrarily generate a string that looks like "p1, p2, p3" for arbitrary N, and parse it into a list of parameters:

params_string = ", ".join("p{}".format(i + 1) for i in range(n))
params = parameters(params_string)

The notion of parsing a string to get a list of parameters sounds smelly to me, and I'm sure there's a nicer way to programatically declare a bunch of parameters, but this will be as close as possible to what your original code was doing.

EDIT: Looking at the Symfit documentation, it seems like parameters(s) is just a shortcut, and you can actually just do:

params = [Parameter("p{}".format(i + 1)) for i in range(n)]

which doesn't require you to build your own joined string of all the parameter names just so Symfit can split them back into individual parameter names. This would also allow you to define other properties for your parameters, such as their initial values or their min/max bounds.

• You can generalize the Eq constraint:

coeffs = [xdata[0]] + [(xdata[i+1] - xdata[i]) for i in range(n-1)]
eq_constraint = Eq(sum(param * coeff for param, coeff in zip(params, coeffs), 1)

Or, as another answer does this, using numpy operations:

coeffs = np.concat(xdata[:1], np.diff(xdata))
eq_constraint = Eq(np.sum(params * coeffs), 1)

• You can generalize the Ge constraints:

ge_constraints = [Ge(params[i + 1], params[i]) for i in range(n - 1)] + [Ge(params[-1], 0]

constraints = [eq_constraint] + ge_constraints

Again this could be done using numpy operations, but I'll leave that to @user3483203's answer.

• You can multiply all the parameters using reduce:

model = reduce(lambda x, y: x * y, params, 1)

Or using numpy.prod:

model = np.prod(params)

This should be enough to generalize the above to arbitrary N.