Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def f2(z, *params):
x, y = z
a, b, c, d, e, f, g, h, i, j, k, l, scale = params
return (-g*np.exp(-((x-h)**2 + (y-i)**2) / scale))
def f3(z, *params):
x, y = z
a, b, c, d, e, f, g, h, i, j, k, l, scale = params
return (-j*np.exp(-((x-k)**2 + (y-l)**2) / scale))
def f(z, *params):
return f1(z, *params) + f2(z, *params) + f3(z, *params)
# setup for scipy-brute optimization #
# setup for lmfit-brute optimization #
params_lmfit = lmfit.Parameters()
params_lmfit.add_many(
('a', 2, False, None, None, None),
('b', 3, False, None, None, None),
('c', 7, False, None, None, None),
('d', 8, False, None, None, None),
('e', 9, False, None, None, None),
('f', 10, False, None, None, None),
('g', 44, False, None, None, None),
('h', -1, False, None, None, None),
('i', 2, False, None, None, None),
('j', 26, False, None, None, None),
('k', 1, False, None, None, None),
('l', -2, False, None, None, None),
('scale', 0.5, False, None, None, None),
('x', -4.0, True, -4.0, 4.0, None, None),
('y', -2.0, True, -2.0, 2.0, None, None),
- label: dataset1
type: spectral
megacomplexes: [mc1]
path: ''
irf: irf1
'''
initial_parameter = [123e-3, 234e-4, 567e-5, -0.05, 0.1]
times = np.concatenate([np.arange(-10, -1, 0.1).flatten(),
np.arange(-1, 10, 0.01).flatten(),
np.arange(10, 50, 1.5).flatten(),
np.arange(100, 1000, 15).flatten()])
wavenum = np.arange(12820, 15120, 4.6)
simparams = Parameters()
simparams.add("p1", 101e-3)
simparams.add("p2", 202e-4)
simparams.add("p3", 505e-5)
simparams.add("p4", 0.04)
simparams.add("p5", 0.5)
simparams.pretty_print()
model = parse_yml(fitspec.format(initial_parameter))
fitmodel = KineticSeparableModel(model)
start = time.perf_counter()
data = fitmodel.eval(simparams, *times, **{'dataset':'dataset1',
'noise':True, 'noise_std_dev':0.000001,
'dataset1_x': wavenum,
'amplitudes':[1, 1, 1],
'locations':[14700, 13515, 14180],
data = (5. * np.sin(2 * x - 0.1) * np.exp(-x*x*0.025) +
np.random.normal(size=len(x), scale=0.2) )
# define objective function: returns the array to be minimized
def fcn2min(params, x, data):
""" model decaying sine wave, subtract data"""
amp = params['amp'].value
shift = params['shift'].value
omega = params['omega'].value
decay = params['decay'].value
model = amp * np.sin(x * omega + shift) * np.exp(-x*x*decay)
return model - data
# create a set of Parameters
params = Parameters()
params.add('amp', value= 10, min=0)
params.add('decay', value= 0.1)
params.add('shift', value= 0.0, min=-np.pi/2., max=np.pi/2)
params.add('omega', value= 3.0)
# do fit, here with leastsq model
result = minimize(fcn2min, params, args=(x, data))
# calculate final result
final = data + result.residual
# write error report
report_fit(params)
# try to plot results
np.random.normal(size=x.size, scale=0.2))
# define objective function: returns the array to be minimized
def fcn2min(params, x, data):
"""Model a decaying sine wave and subtract data."""
amp = params['amp']
shift = params['shift']
omega = params['omega']
decay = params['decay']
model = amp * np.sin(x*omega + shift) * np.exp(-x*x*decay)
return model - data
# create a set of Parameters
params = Parameters()
params.add('amp', value=10, min=0)
params.add('decay', value=0.1)
params.add('shift', value=0.0, min=-np.pi/2., max=np.pi/2.)
params.add('omega', value=3.0)
# do fit, here with the default leastsq algorithm
minner = Minimizer(fcn2min, params, fcn_args=(x, data))
result = minner.minimize()
# calculate final result
final = data + result.residual
# write error report
report_fit(result)
# try to plot results
y[x < start] = smallBlind
y[x > maxEquity] = 0
return y
def plot():
import pylab
pylab.plot(xf, yf, 'ko')
pylab.plot(x, y, 'r')
pylab.show()
xf = [minEquity, max_X]
yf = [bigBlind, maxValue]
self.x = x
# create a set of Parameters
params = Parameters()
params.add('pw', value=pw, vary=False)
params.add('adj1', value=1)
params.add('adj2', value=1)
# do fit, here with leastsq model
result = minimize(fcn2min, params, args=(xf, yf), maxfev=3000)
adj1 = result.params['adj1']
adj2 = result.params['adj2']
# plot results
y = finalFunction(x, adj1, adj2, pw, xf[0], maxEquity, smallBlind)
if pl == True:
plot()
self.y = y
def get_lmfitparm(self):
"""
Generates an lmfit parameter class from the present data set.
The following parameters are used:
self.x : 1d ndarray length N
self.y : 1d ndarray length N
self.fit_weights : 1d ndarray length N
self.fit_bool : 1d ndarray length P, bool
self.fit_parm : 1d ndarray length P, float
"""
params = lmfit.Parameters()
# First, add all fixed parameters
for pp in range(len(self.fit_parm)):
if not self.fit_bool[pp]:
if self.fit_bound[pp][0] == self.fit_bound[pp][1]:
self.fit_bound[pp]=[-np.inf, np.inf]
params.add(lmfit.Parameter(name="parm{:04d}".format(pp),
value=self.fit_parm[pp],
vary=self.fit_bool[pp],
min=self.fit_bound[pp][0],
max=self.fit_bound[pp][1],
)
)
# Second, summarize the constraints in a dictionary, where
# keys are the parameter indexes of varied parameters.
for i in range(3):
nonzero = np.abs(models[i]) > 0
indices_final.append(indices[nonzero])
coefs.append(models[i][nonzero])
indices_final = np.array(indices_final)
coefs = np.array(coefs)
# fit a final overall correction to these based just on a linear model
shapes_pred = np.array([cypoly(params_onehot_save, coef, index)
for coef, index in zip(coefs, indices_final)]).T
after_burners = []
for i in range(3):
y = shapes_save[:, i + 3]
ypred = shapes_pred[:, i]
ysk = model_objs[i][0].predict(Xpoly)
lmparams = lmfit.Parameters()
lmparams.add('m', value=1, vary=True)
lmparams.add('b', value=0, vary=True)
results = lmfit.minimize(linear_resid, lmparams, args=(ypred, y))
logger.info('After burner correction for shape {0}'.format(i))
m, b = results.params.valuesdict().values()
after_burners.append([b, m])
logger.info(lmfit.fit_report(results))
after_burners = np.array(after_burners)
logger.info('Final RMSs: {0}'.format(str(rmss)))
analytic_coefs = {'coefs': coefs, 'indices': indices, 'after_burners': after_burners}
return analytic_coefs
hazard_input_vals = [float(x) for x in hazard_input_vals]
for dx in range(1, len(SYS_DS)):
x_sample = hazard_input_vals
y_sample = pb_exceed[dx]
p0m = np.mean(y_sample)
p0s = np.std(y_sample)
if dx >= 2 and p0m < sys_dmg_fitted_params[dx-1].params['median'].value:
p0m = sys_dmg_fitted_params[dx-1].params['median'].value+0.02
# Fit the dist:
params_est = lmfit.Parameters()
params_est.add('median', value=p0m, min=MIN, max=MAX)
params_est.add('logstd', value=p0s, min=MIN, max=MAX)
params_est.add('loc', value=0.0, vary=False)
sys_dmg_fitted_params[dx] = lmfit.minimize(
res_lognorm_cdf,
params_est,
args=(x_sample, y_sample),
method='leastsq',
nan_policy='omit',
maxfev=1000
)
sys_dmg_model.loc[SYS_DS[dx]] = (
sys_dmg_fitted_params[dx].params['median'].value,
sys_dmg_fitted_params[dx].params['logstd'].value,
sys_dmg_fitted_params[dx].params['loc'].value,
def build_fitmodel(self):
""" use fit components to build model"""
dgroup = self.get_datagroup()
fullmodel = None
params = Parameters()
self.summary = {'components': [], 'options': {}}
for comp in self.fit_components.values():
if comp.usebox is not None and comp.usebox.IsChecked():
for parwids in comp.parwids.values():
params.add(parwids.param)
self.summary['components'].append((comp.mclass.__name__, comp.mclass_kws))
thismodel = comp.mclass(**comp.mclass_kws)
if fullmodel is None:
fullmodel = thismodel
else:
fullmodel += thismodel
self.fit_model = fullmodel
self.fit_params = params
if dgroup is not None:
'min_bound_type', 'max_bound_type')
def observe_changeables(self, changed):
try:
self.lmfit_param.set({changed['name']: changed['value']})
except KeyError:
# the parameter doesn't exist in the lmfit parameter
pass
# print(changed)
class FitModel(Atom):
"""Back-end for the FitController
"""
params = Dict()
lmfit_params = Typed(Parameters)
lmfit_model = Typed(Model)
name = Str()
show_advanced = Bool(False)
show_basic = Bool(True)
def __init__(self, lmfit_model, prefix, name=None):
# print(lmfit_model)
if name is None:
name = str(lmfit_model.name)
self.name = name
if type(lmfit_model) is type:
self.lmfit_model = lmfit_model(name=self.name, prefix=prefix)
else:
self.lmfit_model = lmfit_model
# print(self.lmfit_model.name)
# self.lmfit_model.prefix = "1"