0

I am trying to fit model data (calculated from eR) to my experimental data e_exp. I am not quite sure how to pass constants and variables to func.

import numpy as np
import math
from scipy.optimize import curve_fit, least_squares, minimize

f_exp     = np.array([1, 1.6, 2.7, 4.4, 7.3, 12, 20, 32, 56, 88, 144, 250000])
e_exp     = np.array([7.15, 7.30, 7.20, 7.25, 7.26, 7.28, 7.32, 7.25, 7.35, 7.34, 7.37, 13.55])

ezero     = np.min(e_exp)
einf      = np.max(e_exp)

ig_fc     = 500
ig_alpha  = 0.35

def CCER(einf, ezero, f_exp, fc, alpha):
    x  = [np.log(_ / ig_fc) for _ in f_exp]
    eR = [ezero + 1/2 * (einf - ezero) * (1 + np.sinh((1 - ig_alpha) * _) / (np.cosh((1 - ig_alpha) * _) + np.sin(1/2 * ig_alpha * math.pi))) for _ in x]
    return eR

def func(z):
    return np.sum((CCER(z[0], z[1], z[2], z[3], z[4], z[5]) - e_exp) ** 2)

res = minimize(func, (ig_fc, ig_alpha), method='SLSQP')

einf, ezero, and f_exp are all constant plus the variables I need to optimize are ig_fc and ig_alpha, in which ig stands for initial guess.

How can I make this work?

I am also not sure which of the optimization algorithms from scipy are best suited for my problem (be it curve_fit, least_squares or minimize).

naughty_waves
  • 265
  • 1
  • 13

1 Answers1

1

I believe what you want is the following:

def CCER(x, fc, alpha):
    y = np.log(x/fc)
    eR = ezero + 1/2 * (einf - ezero) * (1 + np.sinh((1 - alpha) * y) / (np.cosh((1 - alpha) * y) + np.sin(1/2 * alpha * math.pi)))
    return eR

res = curve_fit(CCER, f_exp, e_exp, p0=(ig_fc, ig_alpha))

You're passing the first value to CCER as an argument, the two remaining ones (fc and alpha) are then treated as optimizable parameters. All fixed parameters will be read from the outer scope - no need to pass them explicitly to the function here.

Finally, in curve_fit you only need to pass an array of inputs (f_exp) and corresponding outputs (e_exp), as well as - possibly - a tuple of initial guesses p0.

Miłosz Wieczór
  • 131
  • 5
  • Thank you for your answer, @Miłosz Wieczór. I have a bit of trouble understanding exactly what is going on since `einf` and `ezero` are not defined as inputs. Also, are `res[0][0]` and `res[0][1]` the respective `fc` and `alpha` that being subject to optimization? – naughty_waves Jan 19 '20 at 15:32
  • When not found in the local scope of the function, Python will try to search for `einf` and `ezero` in the outer scopes, i.e., the global scope in this case. This is why you don't have to pass them as arguments. Regarding the output, by looking up the documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html you shall see that the first array (`res[0]`) is indeed the list of optimized parameters, and the second (`res[1]`) reports the uncertainties. – Miłosz Wieczór Jan 19 '20 at 20:22
  • Thank you for the explanation, @Miłosz Wieczór! It makes sense to me now. – naughty_waves Jan 20 '20 at 14:44
  • I have a follow-up question for you, @Miłosz Wieczór. What if instead of a vector, `e_exp` is, for example, a 2D array with shape `(2, 12)`? This would yield `einf` and `ezero` not as single values anymore but with the shape `(2,1)`. If I try to run this, I get `error: Result from function call is not a proper array of floats.` This has to have something to do with how Python read the global variables. Is there a way to make it work with multiple `e_exp` so that the calculations are performed on each of them? Thank you, and I apologize for ambushing you like this in the comment section. – naughty_waves Apr 11 '20 at 10:50