-
Notifications
You must be signed in to change notification settings - Fork 95
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
CPU (scipy.optimize.curve fit) and GPU fit results are different. #90
Comments
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. |
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)
|
Don't have time to run this example at the moment, but did you try adding noise to the data? |
I tested this code a long time ago, however, I think the outcome was more inaccurate in the presence of noisy data. |
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. |
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)
The text was updated successfully, but these errors were encountered: