CPU (scipy.optimize.curve fit) and GPU fit results are different.

[po60nani]

Hi,

I compared the results of 2D Gaussian fit in "scipy.optimize.curve fit" with Gpufit for "GAUSS 2D ELLIPTIC." However, as you can see, the outcomes are very different. Would you please let me know if there is any specific configuration that I missed?

GPU-->p0 true 10.00 mean 5.80 std 0.00
GPU-->p1 true 16.00 mean 16.07 std 0.00
GPU-->p2 true 16.00 mean 15.84 std 0.00
GPU-->p3 true 1.80 mean 1.41 std 0.00
GPU-->p4 true 1.50 mean 1.36 std 0.00
GPU-->p5 true 10.00 mean 10.15 std 0.00
CPU-->p0 true 10.00, scipy 10.00
CPU-->p1 true 16.00, scipy 16.00
CPU-->p2 true 16.00, scipy 16.00
CPU-->p3 true 1.80, scipy 1.80
CPU-->p4 true 1.50, scipy 1.50
CPU-->p5 true 10.00, scipy 10.00

scipy:
fit_params, cov_mat = optimize.curve_fit(gaussian_2d, xy_mesh, data, p0=initial_parameters_, maxfev=5000, method='lm')

Gpufit:
parameters, states, chi_squares, number_iterations, execution_time = gf.fit(data, None, model_id, initial_parameters, None, 5000, None, estimator_id, None)

[jkfindeisen]

There is no specific configuration needed in general. Gpufit should be comparable to scipy in accuracy. If you can create a minimal example that I can run, I could look into it a bit more.

[po60nani]

Thank you for responding so quickly. Here's a simple example:

# true parameters
true_parameters = np.array((10, 16, 16, 1.8, 1.5, 10), dtype=np.float32)
initial_parameters_ = np.array((10, 16, 16, 1.3, 1.3, 10), dtype=np.float32)


# generate x and y values
size_x = 32
g = np.arange(size_x)
yi, xi = np.meshgrid(g, g, indexing='ij')
xi = xi.astype(np.float32)
yi = yi.astype(np.float32)

def generate_gauss_2d(p, xi, yi):
    """
    Generates a 2D Gaussian peak.
    http://gpufit.readthedocs.io/en/latest/api.html#gauss-2d

    :param p: Parameters (amplitude, x,y center position, width, width, offset)
    :param xi: x positions
    :param yi: y positions
    :return: The Gaussian 2D peak.
    """

    y = p[5] + p[0] * np.exp(-(((xi - p[1]) ** 2 / ((p[4] ** 2))) + ((yi - p[2]) ** 2 / ((p[3] ** 2)))))
    return y


def gaussian_2d(xy_mesh, amp, xc, yc, sigma_x, sigma_y, b):
    (x, y) = xy_mesh

    # make the 2D Gaussian matrix
    gauss = b + amp * np.exp(-(((x - xc) ** 2 / ((sigma_x ** 2))) + ((y - yc) ** 2 / ((sigma_y ** 2)))))

    # flatten the 2D Gaussian down to 1D
    return np.ravel(gauss)

# generate data
number_points = -1
number_fits = 1
data = generate_gauss_2d(true_parameters, xi, yi)

data = np.reshape(data, (1, number_points))

data = data.astype(np.float32)

t0 = time.time()
fit_params, cov_mat = optimize.curve_fit(gaussian_2d, [yi, xi], np.ravel(data), p0=initial_parameters_, maxfev=5000, method='lm')
t1 = time.time() - t0
data = np.tile(data, (number_fits, 1))

# estimator ID
estimator_id = gf.EstimatorID.MLE

# model ID
model_id = gf.ModelID.GAUSS_2D_ELLIPTIC

# run Gpufit
number_parameters = 6
initial_parameters = np.expand_dims(initial_parameters_, axis=0)
initial_parameters = np.tile(initial_parameters, (number_fits, 1))

parameters, states, chi_squares, number_iterations, execution_time = gf.fit(data, None, model_id, initial_parameters,
                                                                            None, 5000, None, estimator_id, None)

[superchromix]

Don't have time to run this example at the moment, but did you try adding noise to the data?

[po60nani]

I tested this code a long time ago, however, I think the outcome was more inaccurate in the presence of noisy data.

[superchromix]

It's strange, because the LM algorithm implemented in Gpufit has been tested extensively, and is as accurate as scipy, etc. Sometimes optimization routines can fail when there is zero noise present in the data, hence my question.