Multidimensional fit using curve_fit where function is over a grid

Question

I am trying to fit a function that is defined on a grid of x,y. The problem is that curve_fit only allows the data points to be a 1D array, while functions defined on a grid are 2D arrays. The following code illustrates what I mean:

import numpy as np
from scipy.optimize import curve_fit

def myfunc(VARS,a,b):
    X,Y = VARS
    return a*(X**b + Y**b)

x = np.linspace(0,1,100)
y = np.linspace(0,1,100)
z1 = x**2+y**2
z2 = np.zeros([100,100])
for i in range(100):
    for j in range(100):
        z2[i][j] = x[i]**2 + y[j]**2
params1,pcov1 = curve_fit(myfunc, (x,y),z1,p0=(1.,1.))
print "params1: {}".format(params1)
params2,pcov2 = curve_fit(myfunc, (x,y),z2,p0=(1.,1.))

and the output is

params1: [1. 2.]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
ValueError: object too deep for desired array

---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
<ipython-input-867-ac34dcacef8e> in <module>()
     14 params1,pcov1 = curve_fit(myfunc, (x,y),z1,p0=(1.,1.))
     15 print "params1: {}".format(params1)
---> 16 params2,pcov2 = curve_fit(myfunc, (x,y),z2,p0=(1.,1.))

c:\python27\lib\site-packages\scipy\optimize\minpack.pyc in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    740         # Remove full_output from kwargs, otherwise we're passing it in twice.
    741         return_full = kwargs.pop('full_output', False)
--> 742         res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    743         popt, pcov, infodict, errmsg, ier = res
    744         cost = np.sum(infodict['fvec'] ** 2)

c:\python27\lib\site-packages\scipy\optimize\minpack.pyc in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    385             maxfev = 200*(n + 1)
    386         retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
--> 387                                  gtol, maxfev, epsfcn, factor, diag)
    388     else:
    389         if col_deriv:

error: Result from function call is not a proper array of floats.

Ideally, a 2D function is defined over all combinations of x,y, while the 2D function permitted by curve_fit are equivalent to

for i in range(100):
    z2[i][i] = x[i]**2 + y[i]**2

so the function is not defined for for example x=0, y=0.5.

Is there a way to fit on a grid?

Answer 1

I have found a solution. I created copies of the values of x using repeat, and then created hard copies of y that correspond to this.

import numpy as np
from scipy.optimize import curve_fit

def myfunc(VARS,a,b):
    X,Y = VARS

    return a*(X**b + Y**b)

x = np.linspace(0,1,100)
y = np.linspace(0,1,100)
x_rep = np.repeat(x,100)
y_rep = np.array([y for i in range(100)]).flatten()


z1 = x**2+y**2
z2 = x_rep**2 + y_rep**2

params1,pcov1 = curve_fit(myfunc, (x,y),z1,p0=(1.,1.))
print "params1: {}".format(params1)
params2,pcov2 = curve_fit(myfunc, (x_rep,y_rep),z2,p0=(1.,1.))
print "params2: {}".format(params2)

Output:
params1: [1. 2.]
params2: [1. 2.]