Too many parameters?

[GORosin]

I'm sorry if this is a stupid question, but I'm getting an error:

  File "/anaconda3/lib/python3.6/collections/__init__.py", line 434, in namedtuple
    exec(class_definition, namespace)
  File "<string>", line 12
SyntaxError: more than 255 arguments

My model looks like this: model_dict={ y1:jim_model(x1,norm11,nthr11,holes_k11,holes_tau11,interface_k11,interface_tau), y2:jim_model(x1,norm12,nthr12,holes_k12,holes_tau12,interface_k12,interface_tau) } fit = Fit(model_dict,x1=x1_data,y1=y1_data,y2=y2_data) fit_results=fit.execute()

where jim_model is a function I defined earlier. The fit works well with just one of the data sets but gives me the above error with two datasets simultaneously. Is this a case of too many parameters at once and might there be a workaround?

[tBuLi]

Interested by problem, can you provide the full traceback so I can see which lines cause problems? Also, is your jim_model a symbolic expression or another Model instance?

[GORosin]

Thanks for replying so quickly :) Heres the full traceback

Traceback (most recent call last):
  File "radiation_fits.py", line 221, in <module>
    plot_and_fit(end_time=end)
  File "radiation_fits.py", line 86, in plot_and_fit
    total_fit(*x_tuple,*y_tuple)
  File "radiation_fits.py", line 163, in total_fit
    fit_results=fit.execute()
  File "/anaconda3/lib/python3.6/site-packages/symfit/core/fit.py", line 581, in execute
    dict(zip(self.model.params, minimizer_ans._popt))
  File "/anaconda3/lib/python3.6/site-packages/symfit/core/fit.py", line 280, in covariance_matrix
    objective=self.objective)
  File "/anaconda3/lib/python3.6/site-packages/symfit/core/fit.py", line 237, in _covariance_matrix
    hess = objective.eval_hessian(**key2str(best_fit_params))
  File "/anaconda3/lib/python3.6/site-packages/symfit/core/objectives.py", line 370, in eval_hessian
    ordered_parameters, **parameters
  File "/anaconda3/lib/python3.6/site-packages/symfit/core/objectives.py", line 223, in eval_hessian
    result = self.model.eval_hessian(**key2str(parameters))._asdict()
  File "/anaconda3/lib/python3.6/site-packages/symfit/core/models.py", line 871, in eval_hessian
    eval_hess_dict = self.hessian_model(*args, **kwargs)._asdict()
  File "/anaconda3/lib/python3.6/site-packages/symfit/core/models.py", line 631, in __call__
    Ans = variabletuple('Ans', self)
  File "/anaconda3/lib/python3.6/site-packages/symfit/core/models.py", line 51, in variabletuple
    named = namedtuple(typename, field_names, *args, **kwargs)
  File "/anaconda3/lib/python3.6/collections/__init__.py", line 434, in namedtuple
    exec(class_definition, namespace)
  File "<string>", line 12
SyntaxError: more than 255 arguments

jim_model is defined seperatley earlier in the script (Jim is our theorist)

def jim_model(x,norm,n_thr,holes_k,holes_tau,interface_k,interface_tau):
    positive_states=holes_k*(1-exp(-x/holes_tau))
    negative_states=interface_k*(1-exp(-x/interface_tau))
    net_charge=positive_states-negative_states-n_thr
    current= norm*(net_charge)**2
    leaked_current = Piecewise((current,net_charge > 0),(0,net_charge<=0))
    return leaked_current

Again, this worked really well for each of the datasets individually.

[tBuLi]

Interesting issue, can you also share how you defined all your parameters and variables? I'm very surprised that it says you have more than 255 since your function definition is still not close to that.

If you have access to a debugger, could you set a breakpoint at lines like

File "/anaconda3/lib/python3.6/site-packages/symfit/core/objectives.py", line 223, in eval_hessian
    result = self.model.eval_hessian(**key2str(parameters))._asdict()

and check what parameters are?

[GORosin]

I don't think there's anything strange when I'm defining the variables

    x1=Variable('x1')
    y1=Variable('y1')
    x2=Variable('x2')
    y2=Variable('y2')
    #common parameter(s)
    interface_tau = Parameter('interface_tau')

    #unique parameters
    norm11=Parameter("norm11")
    nthr11=Parameter("nthr11")
    holes_k11=Parameter("holes_k11")
    holes_tau11=Parameter("holes_tau11")
    interface_k11=Parameter("interface_k11")

    norm12=Parameter("norm12")
    nthr12=Parameter("nthr12")
    holes_k12=Parameter("holes_k12")
    holes_tau12=Parameter("holes_tau12")
    interface_k12=Parameter("interface_k12")


    model_dict={y1:jim_model_sym(x1,norm11,nthr11,holes_k11,holes_tau11,interface_k11,interface_tau),
              y2:jim_model_sym(x1,norm12,nthr12,holes_k12,holes_tau12,interface_k12,interface_tau)
    }
    fit = Fit(model_dict,x1=x1_data[0:7000],y1=y1_data[:7000],y2=y2_data[0:7000])
    fit_results=fit.execute()

The parameters at objectives.py are

{'holes_k11': 1.0, 'holes_k12': 1.0, 'holes_tau11': 1.0, 'holes_tau12': 1.0, 'interface_k11': 1.0, 'interface_k12': 1.0, 'interface_tau': 1.0, 'norm11': 1.0, 'norm12': 1.0, 'nthr11': 1.0, 'nthr12': 1.0, 'x1': array([5.21000e-01, 9.98000e-01, 1.47600e+00, ..., 8.32013e+02,
     8.32038e+02, 8.32062e+02])}

I think the problem is in line 434(ish) of __init__.py where the class template it creates in this line

 class_definition = _class_template.format(
        typename = typename,
        field_names = tuple(field_names),
        num_fields = len(field_names),
        arg_list = repr(tuple(field_names)).replace("'", "")[1:-1],
        repr_fmt = ', '.join(_repr_template.format(name=name)
                             for name in field_names),
        field_defs = '\n'.join(_field_template.format(index=index, name=name)
                               for index, name in enumerate(field_names))
    )

The __new__() method of this class looks lilke this

def __new__(_cls,  d2y1_dholes_k112, d2y1_dholes_k11dholes_k12, 
d2y1_dholes_k11dholes_tau11, d2y1_dholes_k11dholes_tau12, d2y1_dholes_k11dinterface_k11, 
d2y1_dholes_k11dinterface_k12, d2y1_dholes_k11dinterface_tau, d2y1_dholes_k11dnorm11, 
d2y1_dholes_k11dnorm12, d2y1_dholes_k11dnthr11, d2y1_dholes_k11dnthr12, d2y1_dholes_k122, 
d2y1_dholes_k12dholes_k11, d2y1_dholes_k12dholes_tau11, d2y1_dholes_k12dholes_tau12, 
d2y1_dholes_k12dinterface_k11, d2y1_dholes_k12dinterface_k12, d2y1_dholes_k12dinterface_tau, 
d2y1_dholes_k12dnorm11, d2y1_dholes_k12dnorm12, d2y1_dholes_k12dnthr11, 
d2y1_dholes_k12dnthr12, d2y1_dholes_tau112, d2y1_dholes_tau11dholes_k11, 
d2y1_dholes_tau11dholes_k12, d2y1_dholes_tau11dholes_tau12, 
d2y1_dholes_tau11dinterface_k11, d2y1_dholes_tau11dinterface_k12, 
d2y1_dholes_tau11dinterface_tau, d2y1_dholes_tau11dnorm11, d2y1_dholes_tau11dnorm12, 
d2y1_dholes_tau11dnthr11, d2y1_dholes_tau11dnthr12, d2y1_dholes_tau122, 
....
about 100 lines of this pattern
....
dy2_dholes_k12, dy2_dholes_tau11, dy2_dholes_tau12, dy2_dinterface_k11, dy2_dinterface_k12, 
dy2_dinterface_tau, dy2_dnorm11, dy2_dnorm12, dy2_dnthr11, dy2_dnthr12, y1, y2):
    

Which has more than 255 parameters (313 I think).

[tBuLi]

Sorry for the long silence, I was to shocked to respond. That's a very serious issue, thanks for reporting it! As a short term solution, you can probably use GradientModel or even CallableModel instead of the default Model like so:

model_dict={y1:jim_model_sym(x1,norm11,nthr11,holes_k11,holes_tau11,interface_k11,interface_tau),
              y2:jim_model_sym(x1,norm12,nthr12,holes_k12,holes_tau12,interface_k12,interface_tau)
}
model = GradientModel(model_dict)
fit = Fit(model, x1=x1_data[0:7000],y1=y1_data[:7000],y2=y2_data[0:7000])

This way the Hessian will not be computed, and therefore the problem should not occur. Only downside is that now you can't use Hessian methods, but since we currently only have one of those that is probably not an issue.

I'll need to think about how to solve this issue properly but at least Jim should be happy.

[GORosin]

Thank you very much, I got that to work eventually. Interestingly, the issue is actually solved in python 3.7, because they did away with the 255 parameter limit. For python 3.6 and below, I think it mainly comes down to namedtuple having 255 arguments limit, maybe its possible to replace namedtuple with a dictionary?

Anyway, on my end everything works perfectly after switching to 3.7 so thanks a lot for this awesome project!

[tBuLi]

A dictionary might be a solution, but ultimately I think the solution is the introduction of indexed variables into symfit. With that you could instead write:

    model_dict={y[i]: jim_model_sym(x[i],norm1[i],nthr1[i],holes_k1[i],holes_tau1[i],interface_k1[i],interface_tau)
    }

see here for more details: https://github.com/tBuLi/symfit/pull/179

However, this remains a complicated issue on which I haven't reached consensus with myself on how best to implement it ;).

[lesshaste]

I have the same problem. In my case the large number of arguments seem to be created internally by symfit. Here is my full code:

from symfit import Parameter, Variable, Fit, Eq, Ge, exp, parameters
from symfit.core.objectives import LogLikelihood
import scipy.stats

x = Variable('x')
N = 8
lambdas = parameters(', '.join(['lambdas{}'.format(idx) for idx in range(N)]))
ps = parameters(', '.join(['ps{}'.format(idx) for idx in range(N)]))

model = 0
for idx in range(N):
    ps[idx].min = 0.00001
    lambdas[idx].min = 0.01
    model += ps[idx] * lambdas[idx]*exp(-lambdas[idx]*x)
    
xs = scipy.stats.weibull_min.rvs(loc=0, c=0.5, size=50000)
constraints = [Eq(sum(ps), 1)]
fit = Fit(model, xs, constraints=constraints, objective=LogLikelihood)
fit_result = fit.execute()

This gives:

SyntaxError: more than 255 arguments

I am using Python 3.5

[pckroon]

With the merge of #273 this should be (mostly) fixed (at least, the symptoms alleviated) in master. You could try that. You will have add a minimizer=NelderMead (or Powell) to your Fit call.

[lesshaste]

Unfortunately neither NelderMead nor Powell work well if I use them for LogLikelihood maximisation. Using the code above but setting np.random.seed(333), N = 7 and xs = scipy.stats.weibull_min.rvs(loc=0, c=0.5, size=20000), this is what I get:

The default optimizer runs in 6 seconds and reports no warnings. However:

NelderMead:

The optimizer runs in 20 seconds and outputs:

:2: RuntimeWarning: overflow encountered in multiply :2: RuntimeWarning: overflow encountered in exp :2: RuntimeWarning: invalid value encountered in add :2: RuntimeWarning: overflow encountered in multiply :2: RuntimeWarning: overflow encountered in multiply :2: RuntimeWarning: overflow encountered in multiply :2: RuntimeWarning: overflow encountered in multiply :2: RuntimeWarning: overflow encountered in multiply :2: RuntimeWarning: overflow encountered in multiply :2: RuntimeWarning: overflow encountered in multiply :2: RuntimeWarning: overflow encountered in multiply

Powell:

The optimizer runs in 11 seconds and outputs

/home/user/.local/lib/python3.5/site-packages/scipy/optimize/optimize.py:2340: RuntimeWarning: overflow encountered in double_scalars denom = 2.0 * val :2: RuntimeWarning: overflow encountered in double_scalars :2: RuntimeWarning: invalid value encountered in multiply :2: RuntimeWarning: overflow encountered in double_scalars :2: RuntimeWarning: invalid value encountered in multiply :2: RuntimeWarning: overflow encountered in double_scalars :2: RuntimeWarning: invalid value encountered in multiply :2: RuntimeWarning: overflow encountered in double_scalars :2: RuntimeWarning: invalid value encountered in multiply :2: RuntimeWarning: overflow encountered in double_scalars :2: RuntimeWarning: invalid value encountered in multiply :2: RuntimeWarning: overflow encountered in double_scalars :2: RuntimeWarning: invalid value encountered in multiply :2: RuntimeWarning: overflow encountered in multiply

[pckroon]

The good news is then that it runs :) How is the fit they produced? The RuntimeWarnings are something to keep in mind, but not necessarily a major problem. I did read your original problem too quickly, neither NelderMead, nor Powell will actually respect/use your constraint. You can try COBYLA, SLSQP, or TrustConstr. I'm not sure COBYLA supports equality constraints though, but it should complain loudly if it doesn't. I expect you'll run into the original issue again with TrustConstr. I originally pointed you to NelderMead and Powell, since neither of those use the Jacobian, meaning you'll only have k parameters for the minimizer call (with k the number of Parameters you have). COBYLA and SLSQP use the Jacobian, meaning you'll end up with n*k + k (for n dependent variables in your model) parameters in the minimizer call. TrustConstr will also use the Hessian, resulting in n*k**2 + n*k + k parameters, probably triggering your original error.

[lesshaste]

The example I showed was with N reduced to 7 so that it runs with symfit 0.51. I haven't tried symfit master with the new fix included which is what would be needed for N = 8 (that's where the too many parameters error kicks in).

You are right COBYLA doesn't support Eq.

SLSQP works perfectly (for N=7). Should it now work for N=8 with the fix in master?

TrustConstr (which I may not know how to use.. I just blindly put in minimizer=TrustConstr) takes so long it hasn't finished yet.

[pckroon]

SLSQP works perfectly (for N=7). Should it now work for N=8 with the fix in master?

Yes. The following table lists your N against the number of parameters in your model, the number of parameters needed to call your model (assuming 1 component), and the addition of the Jacobian, and the addition of both the Jacobian and the Hessian. It shows nicely the cutoff between N=7 and N=8.

N	N params	model call	+ Jac	+ Jac + Hess
1	2	2	4	8
2	4	4	8	24
3	6	6	12	48
4	8	8	16	80
5	10	10	20	120
6	12	12	24	168
7	14	14	28	224
8	16	16	32	288
9	18	18	36	360
10	20	20	40	440

TrustConstr (which I may not know how to use.. I just blindly put in minimizer=TrustConstr) takes so long it hasn't finished yet.

It should be that easy yes :) It is a more expensive/robust method, so I'm not really surprised it takes longer. I'm not sure on the default tolerances for convergence and constraint satisfaction, so you can try to play with those. All the options specified in [1] you can give to your execute call.

[1] https://docs.scipy.org/doc/scipy/reference/optimize.minimize-trustconstr.html#optimize-minimize-trustconstr

[lesshaste]

That's a very nice table, thank you. It's a shame that with the breadth of optimizers only SLSQP seems to work at all for my simple problem.

Out of interest, is there a simple table setting out which optimizers support which types of constraints and which use the Jacobian, Hessian or neither? I find it hard to dig it out from the scipy docs. For example https://docs.scipy.org/doc/scipy/reference/optimize.minimize-cobyla.html doesn't tell me that COBYLA does not support equality constraints but I found a comment saying it in https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html.

[pckroon]

[2] has a reasonable description of the algorithms. I think it used to mention that COBYLA only does inequality in the past. That said, for the Symfit minimizers you can get the following from our source code (I did reorder them a bit for you):

class NelderMead(ScipyMinimize, BaseMinimizer):
class Powell(ScipyMinimize, BaseMinimizer):

class BFGS(ScipyGradientMinimize):
class LBFGSB(ScipyGradientMinimize, ScipyBoundedMinimizer):

class COBYLA(ScipyConstrainedMinimize, BaseMinimizer):
class SLSQP(ScipyGradientMinimize, ScipyConstrainedMinimize, ScipyBoundedMinimizer):

class TrustConstr(ScipyHessianMinimize, ScipyConstrainedMinimize, ScipyBoundedMinimizer):

# These are special snowflakes
class DifferentialEvolution(ScipyBoundedMinimizer, GlobalMinimizer):
class BasinHopping(ScipyMinimize, GlobalMinimizer):
class MINPACK(ScipyBoundedMinimizer, GradientMinimizer):

So apparently I was wrong before: COBYLA doesn't use any gradient information. If you have constraints Fit will pick SLSQP for you; otherwise it'll take either BFGS or LBFGSB depending on whether your parameters have bounds. Powell and in particular NelderMead we included in Symfit since they're rather robust, although their performance isn't always stellar. TrustConstr I added because Scipy advertises it as their most powerful minimizer, but we probably don't use it to its full potential (since we only use NonlinearConstraints, rather than checking, would be a useful PR, but not fun to do).

You can actually try to use TrustConstr passing hessian='cs' to execute to make it use a finite approximation of the Hessian, instead of the analytical.

[2] https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize

[lesshaste]

Thank you for this.

So if I understand correctly, TrustConstr is in fact the only optimizer that uses the Hessian and this will still fail with the "too many parameters" error after the current fix has been applied (unless we use the hessian='cs' option to stop it using the symbolic Hessian)?

[pckroon]

Yes.

[tBuLi]

If I understand the problem correctly, it comes down to the fact that model calls return a namedtuple, which is limited to 255 in older python versions. I would therefore propose replacing it by a lightweight custom container object instead. I'm not in favor of replacing it by a dict, because that would mess up the current API.

So this new object would just have to allow for both indexing and attribute access just like a namedtuple, but perhaps more importantly it should not be initiated using the **kwargs syntax since this is what causes the issue in the first place. Instead, it only needs to accept two things: variables and output for each of those variables. So I'm thinking of something like this:

class ModelOutput:
    def __init__(self, output_dict):
        """
        :param output_dict: dict of variables -> output.
        """
        if isinstance(output_dict, OrderedDict):
            self.output_dict = output_dict
        else:
            raise TypeError('`input_dict` has to be an OrderedDict.')
        self.variable_names = {var.name: var for var in output_dict}

    def __getattr__(self, name):
        var = self.variable_names[name]
        return self.output_dict[var]
    
    def __getitem__(self, key):
        return self.output_dict.values()[key]

This is an untested draft, but would something like this be a good solution?

[pckroon]

That would work. Alternatively we can stick with the namedtuple, but init it with a list or tuple, rather than **kwargs. That would require careful sorting though, my preference is your ModelOutput. And we need a test for this. PS. Note that you have a .values() too many in __getitem__. In addition, it would be better to do def __getattr__(self, key): return self[key], and make __getitem__ a bit more robust towards input: var = self.variable_names.get(name, name). Or something along those lines.

[tBuLi]

A namedtuple has to be initiated using args and kwargs right, so I don't think that will work.

I'll start coding something up and see about your other suggestions :).

[pckroon]

A namedtuple is a tuple with some extras on the side, so I expect you'll be able to do nt(list).

[tBuLi]

Yes, but then you wont have the nice __getattr__ anymore, then you can only access the full list, not the entries within that list as we need. And inheriting from a namedtuple to solve that is not really possible becuase namedtuple is actually a factory function. So a custom object does seem like the best solution, I'm almost done with making a PR :).

[pckroon]

Bah, you're right. I'm disappointed about namedtuple now. It seems you can't even really subclass them.