Loading peak_integration.py +93 −58 Original line number Diff line number Diff line Loading @@ -1262,6 +1262,83 @@ def numerical_hessian(x, fun, *args, **kwargs): ############################################################################### ############################################################################### # Loss functions class Loss(object): def __init__(self, bkgr_fun, peak_fun, points, data): self.func_calls = 0 self.grad_calls = 0 self.hess_calls = 0 # background and peak functions self.bkgr_fun = bkgr_fun self.peak_fun = peak_fun # number of background and peak parameters self.nbkgr = self.bkgr_fun.nparams self.npeak = self.peak_fun.nparams # points and data to fit self.points = points self.data = data def fit(self, params): self._fit = self.bkgr_fun(params[:self.nbkgr], self.points) + self.peak_fun(params[self.nbkgr:], self.points) return self._fit def __call__(self, params, *args, **kwargs): self.func_calls+=1 self.func_params = np.asarray(params).copy() return self.loss(params) def gradient(self, params, *args, **kwargs): self.grad_calls += 1 self.grad_params = np.asarray(params).copy() if not hasattr(self, 'func_params') or not np.all(self.grad_params==self.func_params): g = self.__call__(params) dloss_dbkgr = self.bkgr_fun.gradient(params[:self.nbkgr], self.points, dloss_dfit=self.dloss_dfit) dloss_dpeak = self.peak_fun.gradient(params[self.nbkgr:], self.points, dloss_dfit=self.dloss_dfit) return np.concatenate((dloss_dbkgr,dloss_dpeak)) def hessian(self, params, *args, **kwargs): self.hess_calls += 1 self.hess_params = np.asarray(params).copy() if not hasattr(self, 'func_params') or not np.all(self.hess_params==self.func_params): g = self.__call__(params) if not hasattr(self, 'grad_params') or not np.all(self.hess_params==self.grad_params): dg = self.gradient(params) d2loss_dfit2 = self.d2loss_dfit2 d2loss_dbkgr2 = self.bkgr_fun.hessian(params[:self.nbkgr], self.points, dloss_dfit=self._dloss_dfit, d2loss_dfit2=d2loss_dfit2) d2loss_dpeak2 = self.peak_fun.hessian(params[self.nbkgr:], self.points, dloss_dfit=self._dloss_dfit, d2loss_dfit2=d2loss_dfit2) d2loss_dbkgr_dpeak = ((d2loss_dfit2.reshape((-1,1,1)) * self.bkgr_fun.grad[:,:,np.newaxis]) * self.peak_fun.grad[:,np.newaxis,:]).sum(axis=0) d2loss = np.block([ [d2loss_dbkgr2, d2loss_dbkgr_dpeak], [d2loss_dbkgr_dpeak.T, d2loss_dpeak2 ] ]) return d2loss class MLELoss(Loss): def loss(self, params): fit = self.fit(params) return (fit-self.data*np.log(fit)).sum() @property def dloss_dfit(self): self._dloss_dfit = 1 - self.data / self._fit return self._dloss_dfit @property def d2loss_dfit2(self): return self.data / self._fit**2 ############################################################################### ############################################################################### class PeakHistogram(object): Loading Loading @@ -1841,74 +1918,32 @@ class PeakHistogram(object): result = least_squares(residual, #jac=jacobian_residual, x0=params_init, bounds=[params_lbnd,params_ubnd], method='trf', verbose=0, max_nfev=1000) elif loss=='mle': class MLELoss(object): def __init__(self, bkgr_fun, peak_fun): self.func_calls = 0 self.grad_calls = 0 self.hess_calls = 0 # background and peak functions self.bkgr_fun = bkgr_fun self.peak_fun = peak_fun # number of background and peak parameters self.nbkgr = self.bkgr_fun.nparams self.npeak = self.peak_fun.nparams def __call__(self, params): self.func_calls+=1 self.func_params = np.asarray(params).copy() # cache fit for reuse in derivative evaluations self.fit = self.bkgr_fun(params[:self.nbkgr], fit_points) + self.peak_fun(params[self.nbkgr:], fit_points) return (self.fit-fit_data*np.log(self.fit)).sum() def gradient(self, params, *args, **kwargs): self.grad_calls += 1 self.grad_params = np.asarray(params).copy() if not hasattr(self, 'func_params') or not np.all(self.grad_params==self.func_params): g = self.__call__(params) self.dloss_dfit = 1 - fit_data/self.fit dloss_dbkgr = self.bkgr_fun.gradient(params[:self.nbkgr], fit_points, dloss_dfit=self.dloss_dfit) dloss_dpeak = self.peak_fun.gradient(params[self.nbkgr:], fit_points, dloss_dfit=self.dloss_dfit) return np.concatenate((dloss_dbkgr,dloss_dpeak)) def hessian(self, params, *args, **kwargs): self.hess_calls += 1 self.hess_params = np.asarray(params).copy() if not hasattr(self, 'func_params') or not np.all(self.hess_params==self.func_params): g = self.__call__(params) if not hasattr(self, 'grad_params') or not np.all(self.hess_params==self.grad_params): dg = self.gradient(params) d2loss_dfit2 = fit_data / self.fit**2 d2loss_dbkgr2 = self.bkgr_fun.hessian(params[:self.nbkgr], fit_points, dloss_dfit=self.dloss_dfit, d2loss_dfit2=d2loss_dfit2) d2loss_dpeak2 = self.peak_fun.hessian(params[self.nbkgr:], fit_points, dloss_dfit=self.dloss_dfit, d2loss_dfit2=d2loss_dfit2) d2loss_dbkgr_dpeak = ((d2loss_dfit2.reshape((-1,1,1)) * self.bkgr_fun.grad[:,:,np.newaxis]) * self.peak_fun.grad[:,np.newaxis,:]).sum(axis=0) d2loss = np.block([ [d2loss_dbkgr2, d2loss_dbkgr_dpeak], [d2loss_dbkgr_dpeak.T, d2loss_dpeak2 ] ]) return d2loss loss_fun = MLELoss(bkgr_fun=self.bkgr_fun, peak_fun=self.peak_fun, points=fit_points, data=fit_data) # from scipy.optimize import BFGS # class myBFGS(BFGS): # def initialize(self, n, approx_type): # super().initialize(n, approx_type) # if self.approx_type == 'hess': # self.B = MLELoss().hessian(params_init) # # self.B = np.eye(n, dtype=float) # self.B = loss_fun.hessian(params_init) # else: # self.H = np.linalg.inv(MLELoss().hessian(params_init)) # # self.H = np.eye(n, dtype=float) loss_fun = MLELoss(bkgr_fun=self.bkgr_fun, peak_fun=self.peak_fun) # self.H = np.linalg.inv(loss_fun.hessian(params_init)) # result = fmin_bfgs_2(loss_fun, params_init, fprime=loss_fun.gradient, gtol=1e-6, maxiter=1000, disp=2, init_hess=None)#(loss_fun.hessian(params_init)+loss_fun.hessian(params_init).T)/2) # result = fmin_bfgs_2(loss_fun, result.x, fprime=loss_fun.gradient, gtol=1e-6, maxiter=1000, disp=2, init_hess=loss_fun.hessian(result.x)) # result = minimize(loss_fun, # jac = loss_fun.gradient, # hess = loss_fun.hessian, # x0=result.x, bounds=tuple(zip(params_lbnd,params_ubnd)), # method='Newton-CG', options={'maxiter':1000, 'maxfun':1000, 'maxls':20, 'maxcor':100, 'ftol':1.e-8, 'gtol':1.e-6, 'xtol':1.e-6, 'disp':True} #_debug} # ) result = minimize(loss_fun, jac = loss_fun.gradient, hess = None, #loss_fun.hessian, hess = loss_fun.hessian, # hess = myBFGS(), #loss_fun.hessian, x0=params_init, bounds=tuple(zip(params_lbnd,params_ubnd)), method=solver, options={'maxiter':100, 'maxfun':1000, 'maxls':20, 'maxcor':100, 'ftol':1.e-8, 'gtol':1.e-6, 'xtol':1.e-6, 'disp':_debug} method=solver, options={'maxiter':1000, 'maxfun':1000, 'maxls':20, 'maxcor':100, 'ftol':1.e-8, 'gtol':1.e-6, 'xtol':1.e-6, 'disp':True} #_debug} # method='L-BFGS-B', options={'maxiter':1000, 'maxfun':1000, 'maxls':20, 'maxcor':100, 'ftol':1.e-6, 'gtol':1.e-6, 'disp':True} # method='BFGS', options={'maxiter':1000, 'gtol':1.e-6, 'norm':np.inf, 'disp':True} # method='CG', options={'maxiter':1000, 'gtol':1.e-5, 'norm':np.inf, 'disp':True} Loading Loading
peak_integration.py +93 −58 Original line number Diff line number Diff line Loading @@ -1262,6 +1262,83 @@ def numerical_hessian(x, fun, *args, **kwargs): ############################################################################### ############################################################################### # Loss functions class Loss(object): def __init__(self, bkgr_fun, peak_fun, points, data): self.func_calls = 0 self.grad_calls = 0 self.hess_calls = 0 # background and peak functions self.bkgr_fun = bkgr_fun self.peak_fun = peak_fun # number of background and peak parameters self.nbkgr = self.bkgr_fun.nparams self.npeak = self.peak_fun.nparams # points and data to fit self.points = points self.data = data def fit(self, params): self._fit = self.bkgr_fun(params[:self.nbkgr], self.points) + self.peak_fun(params[self.nbkgr:], self.points) return self._fit def __call__(self, params, *args, **kwargs): self.func_calls+=1 self.func_params = np.asarray(params).copy() return self.loss(params) def gradient(self, params, *args, **kwargs): self.grad_calls += 1 self.grad_params = np.asarray(params).copy() if not hasattr(self, 'func_params') or not np.all(self.grad_params==self.func_params): g = self.__call__(params) dloss_dbkgr = self.bkgr_fun.gradient(params[:self.nbkgr], self.points, dloss_dfit=self.dloss_dfit) dloss_dpeak = self.peak_fun.gradient(params[self.nbkgr:], self.points, dloss_dfit=self.dloss_dfit) return np.concatenate((dloss_dbkgr,dloss_dpeak)) def hessian(self, params, *args, **kwargs): self.hess_calls += 1 self.hess_params = np.asarray(params).copy() if not hasattr(self, 'func_params') or not np.all(self.hess_params==self.func_params): g = self.__call__(params) if not hasattr(self, 'grad_params') or not np.all(self.hess_params==self.grad_params): dg = self.gradient(params) d2loss_dfit2 = self.d2loss_dfit2 d2loss_dbkgr2 = self.bkgr_fun.hessian(params[:self.nbkgr], self.points, dloss_dfit=self._dloss_dfit, d2loss_dfit2=d2loss_dfit2) d2loss_dpeak2 = self.peak_fun.hessian(params[self.nbkgr:], self.points, dloss_dfit=self._dloss_dfit, d2loss_dfit2=d2loss_dfit2) d2loss_dbkgr_dpeak = ((d2loss_dfit2.reshape((-1,1,1)) * self.bkgr_fun.grad[:,:,np.newaxis]) * self.peak_fun.grad[:,np.newaxis,:]).sum(axis=0) d2loss = np.block([ [d2loss_dbkgr2, d2loss_dbkgr_dpeak], [d2loss_dbkgr_dpeak.T, d2loss_dpeak2 ] ]) return d2loss class MLELoss(Loss): def loss(self, params): fit = self.fit(params) return (fit-self.data*np.log(fit)).sum() @property def dloss_dfit(self): self._dloss_dfit = 1 - self.data / self._fit return self._dloss_dfit @property def d2loss_dfit2(self): return self.data / self._fit**2 ############################################################################### ############################################################################### class PeakHistogram(object): Loading Loading @@ -1841,74 +1918,32 @@ class PeakHistogram(object): result = least_squares(residual, #jac=jacobian_residual, x0=params_init, bounds=[params_lbnd,params_ubnd], method='trf', verbose=0, max_nfev=1000) elif loss=='mle': class MLELoss(object): def __init__(self, bkgr_fun, peak_fun): self.func_calls = 0 self.grad_calls = 0 self.hess_calls = 0 # background and peak functions self.bkgr_fun = bkgr_fun self.peak_fun = peak_fun # number of background and peak parameters self.nbkgr = self.bkgr_fun.nparams self.npeak = self.peak_fun.nparams def __call__(self, params): self.func_calls+=1 self.func_params = np.asarray(params).copy() # cache fit for reuse in derivative evaluations self.fit = self.bkgr_fun(params[:self.nbkgr], fit_points) + self.peak_fun(params[self.nbkgr:], fit_points) return (self.fit-fit_data*np.log(self.fit)).sum() def gradient(self, params, *args, **kwargs): self.grad_calls += 1 self.grad_params = np.asarray(params).copy() if not hasattr(self, 'func_params') or not np.all(self.grad_params==self.func_params): g = self.__call__(params) self.dloss_dfit = 1 - fit_data/self.fit dloss_dbkgr = self.bkgr_fun.gradient(params[:self.nbkgr], fit_points, dloss_dfit=self.dloss_dfit) dloss_dpeak = self.peak_fun.gradient(params[self.nbkgr:], fit_points, dloss_dfit=self.dloss_dfit) return np.concatenate((dloss_dbkgr,dloss_dpeak)) def hessian(self, params, *args, **kwargs): self.hess_calls += 1 self.hess_params = np.asarray(params).copy() if not hasattr(self, 'func_params') or not np.all(self.hess_params==self.func_params): g = self.__call__(params) if not hasattr(self, 'grad_params') or not np.all(self.hess_params==self.grad_params): dg = self.gradient(params) d2loss_dfit2 = fit_data / self.fit**2 d2loss_dbkgr2 = self.bkgr_fun.hessian(params[:self.nbkgr], fit_points, dloss_dfit=self.dloss_dfit, d2loss_dfit2=d2loss_dfit2) d2loss_dpeak2 = self.peak_fun.hessian(params[self.nbkgr:], fit_points, dloss_dfit=self.dloss_dfit, d2loss_dfit2=d2loss_dfit2) d2loss_dbkgr_dpeak = ((d2loss_dfit2.reshape((-1,1,1)) * self.bkgr_fun.grad[:,:,np.newaxis]) * self.peak_fun.grad[:,np.newaxis,:]).sum(axis=0) d2loss = np.block([ [d2loss_dbkgr2, d2loss_dbkgr_dpeak], [d2loss_dbkgr_dpeak.T, d2loss_dpeak2 ] ]) return d2loss loss_fun = MLELoss(bkgr_fun=self.bkgr_fun, peak_fun=self.peak_fun, points=fit_points, data=fit_data) # from scipy.optimize import BFGS # class myBFGS(BFGS): # def initialize(self, n, approx_type): # super().initialize(n, approx_type) # if self.approx_type == 'hess': # self.B = MLELoss().hessian(params_init) # # self.B = np.eye(n, dtype=float) # self.B = loss_fun.hessian(params_init) # else: # self.H = np.linalg.inv(MLELoss().hessian(params_init)) # # self.H = np.eye(n, dtype=float) loss_fun = MLELoss(bkgr_fun=self.bkgr_fun, peak_fun=self.peak_fun) # self.H = np.linalg.inv(loss_fun.hessian(params_init)) # result = fmin_bfgs_2(loss_fun, params_init, fprime=loss_fun.gradient, gtol=1e-6, maxiter=1000, disp=2, init_hess=None)#(loss_fun.hessian(params_init)+loss_fun.hessian(params_init).T)/2) # result = fmin_bfgs_2(loss_fun, result.x, fprime=loss_fun.gradient, gtol=1e-6, maxiter=1000, disp=2, init_hess=loss_fun.hessian(result.x)) # result = minimize(loss_fun, # jac = loss_fun.gradient, # hess = loss_fun.hessian, # x0=result.x, bounds=tuple(zip(params_lbnd,params_ubnd)), # method='Newton-CG', options={'maxiter':1000, 'maxfun':1000, 'maxls':20, 'maxcor':100, 'ftol':1.e-8, 'gtol':1.e-6, 'xtol':1.e-6, 'disp':True} #_debug} # ) result = minimize(loss_fun, jac = loss_fun.gradient, hess = None, #loss_fun.hessian, hess = loss_fun.hessian, # hess = myBFGS(), #loss_fun.hessian, x0=params_init, bounds=tuple(zip(params_lbnd,params_ubnd)), method=solver, options={'maxiter':100, 'maxfun':1000, 'maxls':20, 'maxcor':100, 'ftol':1.e-8, 'gtol':1.e-6, 'xtol':1.e-6, 'disp':_debug} method=solver, options={'maxiter':1000, 'maxfun':1000, 'maxls':20, 'maxcor':100, 'ftol':1.e-8, 'gtol':1.e-6, 'xtol':1.e-6, 'disp':True} #_debug} # method='L-BFGS-B', options={'maxiter':1000, 'maxfun':1000, 'maxls':20, 'maxcor':100, 'ftol':1.e-6, 'gtol':1.e-6, 'disp':True} # method='BFGS', options={'maxiter':1000, 'gtol':1.e-6, 'norm':np.inf, 'disp':True} # method='CG', options={'maxiter':1000, 'gtol':1.e-5, 'norm':np.inf, 'disp':True} Loading
