Functions with a variable number of parameters

[GoogleCodeExporter]

Originally submitted by iamholger

Hi, Andy directed me here. 

PyMinuit handles multiple parameters easily - but only if the definition of the
function to minimize is fixed. E.g. a hardcoded f(x,y) works perfectly fine, so 
in our code we have to dynamically generate a string which represents a 
functions
definition of N (N can be anything greater than 0) parameters that is later 
used via
the 'exec' statement to actually define the function PyMinuit is considered to 
optimize.

Original issue reported on code.google.com by jpivar...@gmail.com on 9 Jul 2008 at 12:59

[GoogleCodeExporter]

Right, Minuit doesn't minimize functions with a variable number of parameters.  
I'm having trouble thinking 
of what that means, since I'm thinking of minimization as finding an optimal 
point in an N-dimensional 
space, a mental model that breaks down when the number of dimensions is not 
fixed.

It sounds like you're describing a work-around by constructing the function 
dynamically--- but you mean 
that you construct the function before you minimize, right?  During the 
minimization process, the function is 
constant, and the number of dimensions is constant, right?

If so, then what you describe is not a work-around, but an excellent 
demonstration of how Python provides 
solutions not available in C (or Fortran, Minuit's original implementation), 
where the objective function 
needed to be compiled in; dynamically generating the function at run-time was 
unthinkable.

Ah, I think I'm beginning to see what you wanted, and it's something PyMinuit 
doesn't do but could: an 
objective function defined like f(*args) won't work because PyMinuit inspects 
"f" to determine the argument 
list.  This might be improved by allowing the user to specify the argument list 
exactly when a function is 
generic like this.  I'd rather recommend an idiom on the Python side (as 
opposed to the C side), something 
like

def f(*args):
    return whatever

stringargs = ", ".join(argument_list)
m = minuit.Minuit(eval("lambda %s: f(%s)" % (stringargs, stringargs)))

This still does the work through string processing, but possibly less string 
processing.  It feels a little safer 
than defining the whole function "f" in s string with exec().

Original comment by jpivar...@gmail.com on 9 Jul 2008 at 1:16

• Added labels: Type-Enhancement
• Removed labels: Type-Defect

[GoogleCodeExporter]

Hi - What is the appropriate solution when argument_list contains more than 255
variables?  The code below

    args = ','.join( ['p%d' % (x) for x in range(256)] )
    m = minuit2.Minuit2(eval('lambda %s: fcn(%s)' % (args, args)))

returns the error

   Traceback (most recent call last):
       ...
       m = minuit2.Minuit2(eval('lambda %s: fcn(%s)' % (args, args)))
     File "<string>", line 1
   SyntaxError: more than 255 arguments

Thanks

Original comment by acbec...@gmail.com on 7 Apr 2009 at 9:49

[GoogleCodeExporter]

This is probably a limitation in Python's lambda command.  (I haven't 
experimented,
but I suspect that a function defined by "def" would not have this limitation.) 
Maybe you could use "exec" or even "execfile" to build the function with string
operations, then evaluate it to make it real.

But on the other hand...

Can Minuit even minimize functions with more than 256 parameters?  You guys are
troubling my notions of what Minuit is for and what it can do.  It's a general
non-linear minimizer, so the number of function evaluations needed to find the
minimum grows as a steep polynomial with the number of free parameters.  If your
problem is linear or quadratic, perhaps linear programming or an analytic fit 
would
make the most sense.  With a very large number of dimensions, perhaps a Monte 
Carlo
method would be the most useful, or Monte Carlo to find an approximately 
quadratic
regieme, and then analytic, or find the few parameters that matter and ask 
Minuit to
minimize those.

256 parameters, sheesh!

Original comment by jpivar...@gmail.com on 7 Apr 2009 at 10:43

[GoogleCodeExporter]

[deleted comment]

[GoogleCodeExporter]

I've been  trying to write a general model fit wrapper to pyminuit, so that 
only the
model function (which depends on the input params and xdata) has to be defined. 
But
am finding it extremely limiting  given that the pyminuit object can only be
constructed given its reliance on establishing the fcn input parameters by its
trickery with the  fcn func_code attribute.

I've tried implementing the  lambda function trick inside my own class but I 
can't
get it to work - it complains that "global function fcn"  not found; eval seems 
to
miss the point that my fcn is local, not global! So instead I've had to go 
through
hoops to wrap up the fcn inside another class.

It would be so much easier if  pyminuit accepted a list of parameter names, or 
maybe
a fcn object, which say stored the
parameters internally as a dictionary, and  then we just had to provide a 
method to
return the required fit statistic - i.e. pretty much how the c++ version of 
minuit
works. 

Have you ever considered adding such an alternative way of constructing the 
pyminuit
object?

Original comment by abeardm...@gmail.com on 7 Dec 2009 at 11:18

[GoogleCodeExporter]

First answer:

> Have you ever considered adding such an alternative way of constructing the
pyminuit object?

I hadn't considered it until you brought it up, but it could be good if the 
Minuit
constructor would take a Python callable OR an instance of some Fcn class, 
defined in
the Minuit package.  That way, it would have the current convenience and also 
the
flexibility you describe.

I don't have the time to do that now, but if you're willing to make these 
updates to
a copy of the code, I'll test them and post them as a new version.

Second answer:

If you just want to fool PyMinuit to get something that works, can't you make an
object like the following?

class fake_code:
    def __init__(self):
        self.co_argcount = 3
        self.co_varnames = ("one", "two", "three")

class fake_function:
    def __init__(self):
        self.func_code = fake_code()
    def __call__(self, one, two, three):
        return one**2 + two**2 + three**2

m = minuit.Minuit(fake_function())

The fake_function can have whatever other properties you want.  It is a hacky 
bit of
reverse engineering, but it gets the job done.

Original comment by jpivar...@gmail.com on 7 Dec 2009 at 7:25

[GoogleCodeExporter]

Speaking of which, does the following do what you want?

>>> import minuit
>>> class fake_code:
...     def __init__(self, n):
...         self.co_argcount = n
...         self.co_varnames = tuple(map(str, range(n)))
... 
>>> class fake_function:
...     def __init__(self, n):
...         self.n = n
...         self.func_code = fake_code(n)
...     def __call__(self, *args):
...         output = 0.     
...         for arg in args:
...             output += arg**2
...         return output
... 
>>> m = minuit.Minuit(fake_function(12))
>>> m.values
{'11': 0.0, '10': 0.0, '1': 0.0, '0': 0.0, '3': 0.0, '2': 0.0, '5': 0.0, '4': 
0.0,
'7': 0.0, '6': 0.0, '9': 0.0, '8': 0.0}
>>> for key in m.values.keys():
...     m.values[key] = float(key)
... 
>>> m.values
{'11': 11.0, '10': 10.0, '1': 1.0, '0': 0.0, '3': 3.0, '2': 2.0, '5': 5.0, '4': 
4.0,
'7': 7.0, '6': 6.0, '9': 9.0, '8': 8.0}
>>> m.migrad()
>>> m.values
{'11': 9.4557606189482613e-10, '10': 4.0116692012759358e-07, '1':
8.5257134685434721e-11, '0': 0.0, '3': 2.5834134831370648e-10, '2':
1.7180656897153312e-10, '5': 4.3015990769390555e-10, '4': 
3.4491476341713678e-10,
'7': 6.0199845108854788e-10, '6': 5.1673776368943436e-10, '9':
7.7510620144494169e-10, '8': 6.8855232626674479e-10}

Perhaps we could include your wrapper package in the PyMinuit bundle, so that 
users would

import minuit

for the classic interface and

import multiminuit

for the arbitrary number of parameters version.  Since multiminuit.py would be 
a pure
python package, it could have more features.

Everyone should keep in mind, however, that the underlying Minuit library 
probably
can't hande a very large number of parameters.

Original comment by jpivar...@gmail.com on 7 Dec 2009 at 7:31

[GoogleCodeExporter]

I'd like to push back on this idea that "the underlying Minuit library probably
can't handle a very large number of parameters" since this is the second time 
its
been brought up here.  Do you have any evidence for this?  Its been used in the
particule physics community for many years and I have not heard of any such
limitations, either practical or programmatical.

Original comment by acbec...@gmail.com on 7 Dec 2009 at 7:40

[GoogleCodeExporter]

> I'd like to push back on this idea that "the underlying Minuit library 
probably
> can't handle a very large number of parameters" since this is the second time 
its
> been brought up here.  Do you have any evidence for this?

No, I don't have any evidence for where exactly it breaks down.  (I'm expecting
practical limitations, not programmatic ones.  Of course everything has a 
practical
limitation--- Minuit's limit is just lower than other, more specialized 
methods.  I'd
be very surpised if it can handle non-quadratic functions of a few hundred 
parameters
or more.)

I'm just thinking about how it works (run it with printMode on to see): it has 
to
vary each parameter holding all others fixed, which means that the number of 
function
calls scales exponentially in the number of parameters.  That will quickly 
exhaust
the CPU or memory of the computer you're working on.

If a problem has a very large number of parameters, one tries to cast it as a 
linear
problem and solve it by matrix inversion (polynomial-time in the number of
parameters, if I remember right).  Even that gets unwieldy fairly quickly (at a 
few
thousand parameters or so), and then you have to start thinking about sparse 
matrix
ideas, or about finding a solution to the problem without actually inverting the
matrix (e.g the Millepede II algorithm in detector alignment).

The largest number of free parameters that I've ever heard of people using a 
general
non-linear minimization package for is "dozens", in Dalitz fits of 3-body 
decays. 
For more complicated QCD multi-resonance fitting (like glueball searches), 
people use
something called "wavelet analysis"--- I'm not sure if that's a special 
technique,
more specialized than simple function minimization.

Somewhere in the Minuit documentaion there's a stern warning against using it in
situations where performance is an issue.  It's not meant to solve all 
minimization
problems, it's just a convenient tool for the easy cases.

Original comment by jpivar...@gmail.com on 7 Dec 2009 at 8:13

[GoogleCodeExporter]

Using your extremely useful classes described above, I did a quick test to see 
how
minuit performed when using 256 parameters instead of 12.  When minimizing 
something
quadratic it found the correct minimum to machine precision (1e-14), with 
negligible
impact on system memory, and taking ~5 min on a 3GHz computer.  Minos() also 
found
the correct uncertainties.  

When the exponent was larger than 2 it took far longer (45m), still used 
negligible
memory, and found values that were within the correct ones by 0.1% at most 
(though
the typical difference was an order of magnitude smaller).  I'm not sure I 
entirely
understand the minos() error estimates though.

Overall it seems that Minuit2 can reasonably handle hundreds of parameters.

Original comment by acbec...@gmail.com on 8 Dec 2009 at 7:27

[GoogleCodeExporter]

Apologies for not checking back sooner, but I kinda assumed I would receive a 
gmail 
notification that someone had replied which never happened.

Using the func_code trick, I've managed to get something working reasonably 
well. 
(Though I did discover that an objects func_code is stored in different places 
in 
pytho2.5 and 2.6!). This class generates a generic fcn based on an input model 
and 
statistic function which can then be used in the Minuit constructor :

class FCN:
    """
    Calculates the fitting FCN for pyMinuit(2) given the data (xdata & ydata)
    and model (class Model, with a tuple of initial parameters, params),
    using the class StatFunc to calculate the statistic.
    """
    def __init__(self, xdata, ydata, Model, params, StatFunc):
        self.xdata = xdata
        self.ydata = ydata
        self.Model = Model
        self.model = self.Model(*params)
        self.statfunc = StatFunc()
        # set the func_code to fool the Minuit constructor
        self.func_code = 
fake_code(*self.model.__init__.im_func.func_code.co_varnames[1:])

If there's a more elegant way to do this then let me know.

I've attached an example of use of my pyminuit fitting wrapper code 
(testfit2.py 
which uses fitmodel.py). It's set up to work with numpy and pyminuit2 and 
expects to 
find matplotlib so it can plot the data and best fit model.

I have a question for Jim about use of Minuit.minos(). I usually deal with
chi-square or  -2.0 log likelihood fitting statistics. So for example, a 90% 
confidence error is obtained when the fit statistic increases by 2.706. In the 
c++ 
Minuit2 terminology this is usually referred to as the "up" parameter. 

Now in pyminuit, the minos() method refers to the  "nsigma" for the error 
estimate. 
Is this parameter the "up" parameter? In which case I set it to 2.7 for 90% 
errors.
Or being as a 90% confidence area corresponds to 1.645 sigma, should I set it 
to 
1.645 ?

Original comment by abeardm...@gmail.com on 23 Dec 2009 at 1:44

Attachments:

• testfit2.py
• fitmodel.py

[GoogleCodeExporter]

> I have a question for Jim about use of Minuit.minos(). I usually deal with
> chi-square or  -2.0 log likelihood fitting statistics. So for example, a 90% 
> confidence error is obtained when the fit statistic increases by 2.706. In 
the c++ 
> Minuit2 terminology this is usually referred to as the "up" parameter. 

> Now in pyminuit, the minos() method refers to the  "nsigma" for the error 
estimate. 
> Is this parameter the "up" parameter? In which case I set it to 2.7 for 90% 
errors.
> Or being as a 90% confidence area corresponds to 1.645 sigma, should I set it 
to 
> 1.645 ?

migrad() and minos() both have an "up" parameter which behaves in the way you 
expect:
1.0 for chi^2 and 0.5 for log-likelihood.  (Documented on
http://code.google.com/p/pyminuit/wiki/FunctionReference .)  The additional 
"sigmas"
parameter is only in minos() and contour(), and it multiplies "up", such that
sigmas=1.0 up=1.0 for a chi^2 fit gives you 68% confidence level contours 
(between
the error bars on both sides of the best-fit), sigmas=2.0 up=1.0 for a chi^2 fit
gives you the 95% confidence level, etc.

Its purpose is to allow you to calculate contours for several confidence 
levels, like
1, 2, and 3 sigma, without touching the "up" parameter, decoupling the plotting
procedure from the question of whether you're doing a chi^2 or log likelihood 
fit.

Original comment by jpivar...@gmail.com on 23 Dec 2009 at 2:42