Commit db433cc9 authored by Viktor Reshniak's avatar Viktor Reshniak
Browse files

Update peak_integration.py

parent 4af36f5f

peak_integration.py

+271 −72
Original line number Diff line number Diff line
@@ -933,6 +933,9 @@ class Histogram(object): 
		# histogram array

		self.data = hist_ws.getNumEventsArray().copy()



		# from  scipy.ndimage import gaussian_filter

		# self.data = gaussian_filter(self.data, sigma=2)



		# limits of the box along each dimension

		self.limits = [(hd.getMinimum(), hd.getMaximum()) for hd in (hist_ws.getDimension(i) for i in range(self.ndims))]
@@ -994,7 +997,215 @@ class Histogram(object): 
			return data, points





	def fit(self, bins=None, return_bins=False, loss='mle', covariance_parameterization='givens', params_init=None, params_lbnd=None, params_ubnd=None):

	def fit_two_level(self, min_level=4, return_bins=False, loss='mle', solver0='BFGS', solver='L-BFGS-B', covariance_parameterization='givens', plot_intermediate=False):

		# shape of the largest subhistogram with shape as a power of 2

		shape2 = [2**int(np.log2(self.shape[d])) for d in range(self.ndims)]



		# smallest power of 2 among all dimensions

		minpow2 = min([int(np.log2(self.shape[d])) for d in range(self.ndims)])



		solver1 = solver0



		params = params_lbnd = params_ubnd = None

		for p in [min_level]:

			start = time.time()





			# left, middle (with power of 2 shape) and right bins

			binsl = [(self.shape[d]-shape2[d])//2 for d in range(self.ndims)]

			binsm = [split_bins([shape2[d]],2**p,recursive=False) for d in range(self.ndims)]

			binsr = [(self.shape[d]-shape2[d]) - binsl[d] for d in range(self.ndims)]



			# # combine two middle bins

			# binsm = [b[:len(b)//2-1]+[b[len(b)//2-1]+b[len(b)//2+1]]+b[len(b)//2+1:] for b in binsm]



			# bins at the current level

			bins = [ ([bl] if bl>0 else [])+bm+([br] if br>0 else []) for bl,bm,br in zip(binsl,binsm,binsr)]



			# fit histogram at the current level

			# params, sucess

			output = self.fit(bins=bins, return_bins=return_bins, loss=loss, solver=solver1, covariance_parameterization=covariance_parameterization, params_init=params, params_lbnd=params_lbnd, params_ubnd=params_ubnd)

			params, sucess = output[0], output[1]



			# skip level if fit not successful

			if not sucess and p<minpow2:

				params = params_lbnd = params_ubnd = None

				continue

			solver1 = solver



			nangles   = self.ndims*(self.ndims-1)//2

			sqrtbkgr  = params[0]

			sqrtintst = params[1]

			cnt       = params[2:2+self.ndims]

			angles    = params[2+self.ndims:2+self.ndims+nangles]

			invrads   = params[2+self.ndims+nangles:2+2*self.ndims+nangles]

			# skewness  = params[2+2*ndims+nangles:2+3*ndims+nangles]



			#######################################################################

			# refine search bounds



			# bounds for the center

			dcnt = [ 10*res*2**(minpow2-p) for res in self.resolution] # search radius is 10 voxels at the current level

			cnt_lbnd = [c-dc for c,dc in zip(cnt,dcnt)]

			cnt_ubnd = [c+dc for c,dc in zip(cnt,dcnt)]



			# bounds for the precision matrix angles

			if minpow2>min_level:

				phi0 = np.pi/1

				phi1 = np.pi#/8

				dphi = phi0 + (p-min_level)/(minpow2-min_level) * (phi1-phi0)

			else:

				dphi = np.pi



			# sc = 2 + (p-min_level)/(minpow2-min_level) * (1-2)

			# sc = 2

			# print(sc)



			######################################

			# bounds for the precision matrix



			# bounds on the semiaxes of the `peak_std` ellipsoid: standard deviations (square roots of the eigenvalues) of the covariance matrix

			peak_std = 4

			max_rads = [ 5/6*rad/peak_std for rad in self.radiuses]   # largest  'peak_std' radius is 5/6 of the box radius

			min_rads = [ 1/2*res/peak_std for res in self.resolution] # smallest 'peak_std' radius is 1/2 of the smallest bin size

			# prec_lbnd = [-np.pi]*num_angles + [ 1/r for r in max_rads]

			# prec_ubnd = [ np.pi]*num_angles + [ 1/r for r in min_rads]

			# prec_lbnd = [max(-np.pi,phi-dphi) for phi in angles] + [ r/2.0 for r in invrads]

			prec_lbnd = [max(-np.pi,phi-dphi) for phi in angles] + [ 1/r for r in max_rads] #[ max((self.limits[d][1]-self.limits[d][0])/4/(4/3),invrads[d]/2.0) for d in range(self.ndims)]

			prec_ubnd = [min( np.pi,phi+dphi) for phi in angles] + [ 1/r for r in min_rads] #[ 100*r for r in invrads]



			# bounds for all parameters

			# params_lbnd = [np.abs(sqrtbkgr)/2, np.abs(sqrtintst)/2] + cnt_lbnd + prec_lbnd #+ [0 for sk in skewness]

			# params_ubnd = [2*np.abs(sqrtbkgr), 2*np.abs(sqrtintst)] + cnt_ubnd + prec_ubnd #+ [0 for sk in skewness]

			params_lbnd = [-np.inf, -np.inf] + cnt_lbnd + prec_lbnd #+ [0 for sk in skewness]

			params_ubnd = [ np.inf,  np.inf] + cnt_ubnd + prec_ubnd #+ [0 for sk in skewness]



			# params[0] = np.abs(sqrtbkgr)

			# params[1] = np.abs(sqrtintst)



			# params_lbnd = params_lbnd + list(params[len(params_lbnd):])

			# params_ubnd = params_ubnd + list(params[len(params_ubnd):])

			# params_lbnd = params_ubnd = None



			#######################################################################



			print(f"Fitted histogram with {2**p:3d} bins: {time.time()-start:.3f} seconds")



			# if plot_intermediate:

			# 	plot_fit(hist_ws, params, bins, prefix=f"{p}", peak_id=1074, peak_hkl=[2.0,-2.0,-9.0], peak_std=4, bkgr_std=7, detector_mask=None, log=True)



		start = time.time()

		output = self.fit(return_bins=return_bins, loss=loss, solver=solver, covariance_parameterization=covariance_parameterization, params_init=params, params_lbnd=params_lbnd, params_ubnd=params_ubnd)

		print(f"Fitted histogram with origianl bins: {time.time()-start:.3f} seconds")



		return output





	def fit_multilevel(self, min_level=4, return_bins=False, loss='mle', solver0='BFGS', solver='L-BFGS-B', covariance_parameterization='givens', plot_intermediate=False):

		# shape of the largest subhistogram with shape as a power of 2

		shape2 = [2**int(np.log2(self.shape[d])) for d in range(self.ndims)]



		# smallest power of 2 among all dimensions

		minpow2 = min([int(np.log2(self.shape[d])) for d in range(self.ndims)])



		solver1 = solver0



		params = params_lbnd = params_ubnd = None

		for p in range(min_level,minpow2+1):

			start = time.time()





			# left, middle (with power of 2 shape) and right bins

			binsl = [(self.shape[d]-shape2[d])//2 for d in range(self.ndims)]

			binsm = [split_bins([shape2[d]],2**p,recursive=False) for d in range(self.ndims)]

			binsr = [(self.shape[d]-shape2[d]) - binsl[d] for d in range(self.ndims)]



			# # combine two middle bins

			# binsm = [b[:len(b)//2-1]+[b[len(b)//2-1]+b[len(b)//2+1]]+b[len(b)//2+1:] for b in binsm]



			# bins at the current level

			bins = [ ([bl] if bl>0 else [])+bm+([br] if br>0 else []) for bl,bm,br in zip(binsl,binsm,binsr)]



			# fit histogram at the current level

			# params, sucess

			output = self.fit(bins=bins, return_bins=return_bins, loss=loss, solver=solver1, covariance_parameterization=covariance_parameterization, params_init=params, params_lbnd=params_lbnd, params_ubnd=params_ubnd)

			params, sucess = output[0], output[1]



			# skip level if fit not successful

			# if not sucess and p<minpow2:

			# 	params = params_lbnd = params_ubnd = None

			# 	continue

			solver1 = solver



			nangles   = self.ndims*(self.ndims-1)//2

			sqrtbkgr  = params[0]

			sqrtintst = params[1]

			cnt       = params[2:2+self.ndims]

			angles    = params[2+self.ndims:2+self.ndims+nangles]

			invrads   = params[2+self.ndims+nangles:2+2*self.ndims+nangles]

			# skewness  = params[2+2*ndims+nangles:2+3*ndims+nangles]



			#######################################################################

			# refine search bounds



			# bounds for the center

			dcnt = [ 10*res*2**(minpow2-p) for res in self.resolution] # search radius is 10 voxels at the current level

			cnt_lbnd = [c-dc for c,dc in zip(cnt,dcnt)]

			cnt_ubnd = [c+dc for c,dc in zip(cnt,dcnt)]



			# bounds for the precision matrix angles

			if minpow2>min_level:

				phi0 = np.pi/1

				phi1 = np.pi#/8

				dphi = phi0 + (p-min_level)/(minpow2-min_level) * (phi1-phi0)

			else:

				dphi = np.pi



			# sc = 2 + (p-min_level)/(minpow2-min_level) * (1-2)

			# sc = 2

			# print(sc)



			######################################

			# bounds for the precision matrix



			# bounds on the semiaxes of the `peak_std` ellipsoid: standard deviations (square roots of the eigenvalues) of the covariance matrix

			peak_std = 4

			max_rads = [ 5/6*rad/peak_std for rad in self.radiuses]   # largest  'peak_std' radius is 5/6 of the box radius

			min_rads = [ 1/2*res/peak_std for res in self.resolution] # smallest 'peak_std' radius is 1/2 of the smallest bin size

			# prec_lbnd = [-np.pi]*num_angles + [ 1/r for r in max_rads]

			# prec_ubnd = [ np.pi]*num_angles + [ 1/r for r in min_rads]

			# prec_lbnd = [max(-np.pi,phi-dphi) for phi in angles] + [ r/2.0 for r in invrads]

			prec_lbnd = [max(-np.pi,phi-dphi) for phi in angles] + [ 1/r for r in max_rads] #[ max((self.limits[d][1]-self.limits[d][0])/4/(4/3),invrads[d]/2.0) for d in range(self.ndims)]

			prec_ubnd = [min( np.pi,phi+dphi) for phi in angles] + [ 1/r for r in min_rads] #[ 100*r for r in invrads]



			# bounds for all parameters

			# params_lbnd = [np.abs(sqrtbkgr)/2, np.abs(sqrtintst)/2] + cnt_lbnd + prec_lbnd #+ [0 for sk in skewness]

			# params_ubnd = [2*np.abs(sqrtbkgr), 2*np.abs(sqrtintst)] + cnt_ubnd + prec_ubnd #+ [0 for sk in skewness]

			params_lbnd = [-np.inf, -np.inf] + cnt_lbnd + prec_lbnd #+ [0 for sk in skewness]

			params_ubnd = [ np.inf,  np.inf] + cnt_ubnd + prec_ubnd #+ [0 for sk in skewness]



			# params[0] = np.abs(sqrtbkgr)

			# params[1] = np.abs(sqrtintst)



			# params_lbnd = params_lbnd + list(params[len(params_lbnd):])

			# params_ubnd = params_ubnd + list(params[len(params_ubnd):])

			# params_lbnd = params_ubnd = None



			#######################################################################



			print(f"Fitted histogram with {2**p:3d} bins: {time.time()-start:.3f} seconds")



			# if plot_intermediate:

			# 	plot_fit(hist_ws, params, bins, prefix=f"{p}", peak_id=1074, peak_hkl=[2.0,-2.0,-9.0], peak_std=4, bkgr_std=7, detector_mask=None, log=True)



		# start = time.time()

		# output = self.fit(return_bins=return_bins, loss=loss, solver=solver, covariance_parameterization=covariance_parameterization, params_init=params, params_lbnd=params_lbnd, params_ubnd=params_ubnd)

		# print(f"Fitted histogram with origianl bins: {time.time()-start:.3f} seconds")



		return output





	def fit(self, bins=None, return_bins=False, loss='mle', solver='BFGS', covariance_parameterization='givens', params_init=None, params_lbnd=None, params_ubnd=None):

		'''

		Inputs

		------
@@ -1048,9 +1259,9 @@ class Histogram(object): 
		fit_points = fit_points.reshape((-1,self.ndims))



		# detector_mask = None

		if self.detector_mask is None:

			detector_mask = fit_data==fit_data

		else:

		# if self.detector_mask is None:

		# 	detector_mask = fit_data==fit_data

		if self.detector_mask is not None:

			detector_mask = rebin_histogram(self.detector_mask.astype(int), bins)>0

			fit_data   = fit_data.ravel()[detector_mask.ravel()]

			fit_points = fit_points[detector_mask.ravel(),:]
@@ -1093,8 +1304,8 @@ class Histogram(object): 
				max_rads = [ 5/6*rad/peak_std for rad in self.radiuses]   # largest  'peak_std' radius is 5/6 of the box radius

				min_rads = [ 1/2*res/peak_std for res in self.resolution] # smallest 'peak_std' radius is 1/2 of the smallest bin size

				# `num_angles` angles and `ndims` square roots of the eigenvalues of the precision matrix

				prec_init = [     0]*num_angles + [ 1/r for r in ini_rads]

				# prec_init = [     0]*num_angles + std_1d

				prec_init = [     0]*num_angles + [ 1/r for r in ini_rads]

				prec_lbnd = [-np.pi]*num_angles + [ 1/r for r in max_rads]

				prec_ubnd = [ np.pi]*num_angles + [ 1/r for r in min_rads]

			elif covariance_parameterization=='cholesky':
@@ -1109,10 +1320,10 @@ class Histogram(object): 
				prec_lbnd = [-1000]*(self.ndims**2)

				prec_ubnd = [ 1000]*(self.ndims**2)



			# initialization and bounds for the skewness

			skew_init = [0]*self.ndims

			skew_lbnd = [0]*self.ndims

			skew_ubnd = [0]*self.ndims

			# # initialization and bounds for the skewness

			# skew_init = [0]*self.ndims

			# skew_lbnd = [0]*self.ndims

			# skew_ubnd = [0]*self.ndims



		###################################

		# initialization and bounds for all parameters
@@ -1121,11 +1332,16 @@ class Histogram(object): 
		if params_lbnd is None: params_lbnd = bkgr_lbnd + intst_lbnd + cnt_lbnd + prec_lbnd #+ skew_lbnd

		if params_ubnd is None: params_ubnd = bkgr_ubnd + intst_ubnd + cnt_ubnd + prec_ubnd #+ skew_ubnd



		# params_init = bkgr_init + intst_init + list(params_init[2:])

		# params_lbnd = bkgr_lbnd + intst_lbnd + list(params_lbnd[2:])

		# params_ubnd = bkgr_ubnd + intst_ubnd + list(params_ubnd[2:])





		# params_lbnd = [p-1.e-8 for p in params_init]

		# params_ubnd = [p+1.e-8 for p in params_init]



		# number of background and peak parameters

		nbkgr = len(bkgr_init)

		nbkgr = 1 #len(bkgr_init)

		npeak = len(params_init) - nbkgr



		###########################################################################
@@ -1201,52 +1417,22 @@ class Histogram(object): 
					d2loss[:nbkgr,nbkgr:] = 2*params[0] * (d2loss_dfit2.reshape((-1,1)) * gaussian_peak.dg).sum(axis=0)

					d2loss[nbkgr:,:nbkgr] = d2loss[:nbkgr,nbkgr:].T





					return d2loss

			# @MemoizeJacHess

			# def peak_loss(params):

			# 	# background

			# 	bkgr = params[0]**2



			# 	# peak

			# 	g,dg,d2g = gaussian_mixture(params[nbkgr:],fit_points,npeaks=1,covariance_parameterization=covariance_parameterization,return_gradient=True,return_hessian=True)

			# 	# g = gaussian_mixture(params[nbkgr:],fit_points,npeaks=1,covariance_parameterization=covariance_parameterization,return_gradient=False,return_hessian=False)

			# 	print(g.shape,dg[0].shape)

			# 	exit()

			# 	fit = bkgr + g



			# 	# derivative wrt fit

			# 	dloss_dfit   = 1 - fit_data/fit

			# 	d2loss_dfit2 = fit_data/fit**2



			# 	# gradient of the loss

			# 	dloss_dbkgr = [2*params[0]*dloss_dfit.sum()]#,0,0,0]

			# 	# dloss_dpeak = np.concatenate([ (dloss_dfit.reshape((-1,1))*dgi).sum(axis=0) for dgi in dg ])

			# 	# dloss = np.concatenate((dloss_dbkgr,dloss_dpeak))

			# 	# dloss_dpeak = (dloss_dfit.reshape((-1,1))*dg).sum(axis=0)

			# 	dloss_dpeak = np.einsum('k,ki->i',dloss_dfit,dg)

			# 	dloss = np.concatenate((dloss_dbkgr,dloss_dpeak))

			# 	print(dloss)

			# 	exit()



			# 	# hessian of the loss

			# 	d2loss = np.zeros((len(params),len(params)))

			# 	d2loss_dbkgr2 = 2 * dloss_dfit.sum() + (2*params[0])**2 * d2loss_dfit2.sum()

			# 	d2loss_dpeak2 = (dloss_dfit.reshape((-1,1,1))*d2g + (d2loss_dfit2.reshape((-1,1))*dg)[...,np.newaxis]*dg[:,np.newaxis,:]).sum(axis=0)

			# 	d2loss[:nbkgr,:nbkgr] = d2loss_dbkgr2

			# 	d2loss[nbkgr:,nbkgr:] = d2loss_dpeak2

			# from scipy.optimize import BFGS

			# print(BFGS)

			# exit()



			# 	return (fit-fit_data*np.log(fit)).sum(), dloss, d2loss #None, None

			# peak_loss = MemoizeJacHess(peak_loss)

			peak_loss = MLELoss()

			result = minimize(peak_loss,

				jac  = peak_loss.gradient, #gaussian.gradient, #peak_loss.derivative,

				hess = peak_loss.hessian, #peak_loss.hessian,

				jac  = peak_loss.gradient,

				hess = peak_loss.hessian,

				x0=params_init, bounds=tuple(zip(params_lbnd,params_ubnd)),

				# method='L-BFGS-B', options={'maxiter':1000, 'maxfun':1000, 'maxls':20, 'maxcor':100, 'ftol':1.e-5, 'gtol':1.e-5, 'disp':False}

				# method='BFGS', options={'maxiter':1000, 'gtol':1.e-5, 'norm':np.inf, 'disp':True}

				method=solver, options={'maxiter':100, 'maxfun':1000, 'maxls':20, 'maxcor':100, 'ftol':1.e-8, 'gtol':1.e-6, 'xtol':1.e-6, 'disp':False}

				# 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}

				method='Newton-CG', options={'maxiter':1000, 'xtol':1.e-5, 'disp':True}

				# method='Newton-CG', options={'maxiter':1000, 'xtol':1.e-6, 'disp':True}

				# method='Nelder-Mead', options={'maxiter':1000, 'disp':False}

				# method='TNC', options={'scale':None, 'maxfun':1000, 'ftol':1.e-3, 'gtol':1.e-5, 'disp':True}

				# method='dogleg', options={'maxiter':1000, 'tol':1.e-6, 'gtol':1.e-8, 'disp':True}
@@ -1254,7 +1440,12 @@ class Histogram(object): 
				# method='trust-exact', options={'maxiter':1000, 'gtol':1.e-4, 'disp':True}

				# method='trust-constr', options={'maxiter':1000, 'gtol':1.e-4, 'disp':True}

				)

			# restrict parameters

			if solver!='L-BFGS-B':

				for i in range(len(result.x)):

					result.x[i] = max(min(result.x[i],params_ubnd[i]),params_lbnd[i])



			# print(peak_loss.func_calls, peak_loss.grad_calls, peak_loss.hess_calls)

			# g,dg = gaussian_mixture(result.x[nbkgr:],fit_points,npeaks=1,covariance_parameterization=covariance_parameterization,return_gradient=True)

			# fit = result.x[0]**2 + g.reshape(fit_data.shape)

			# res = (fit-fit_data) / fit
@@ -1264,31 +1455,33 @@ class Histogram(object): 
			# print(result.jac)

			# print(grad)



			np.set_printoptions(precision=1, linewidth=200)



			self.fit_params = result.x

			# fit_params[0] = fit_params[0] / np.sqrt(sc)

			# fit_params[1] = fit_params[1] / np.sqrt(sc)

			# return self.fit_params, True, bins

			# print(result)

			print(result.success)

			# print('Chi2: ',chi2(fit_params))



			# print(np.linalg.eig(peak_loss.hessian(self.fit_params))[0])

			# print(np.linalg.eig(numerical_hessian(self.fit_params, lambda y: peak_loss(y)))[0])



			# # start = time.time()

			# # g,dg,d2g = gaussian_mixture(self.fit_params[nbkgr:],fit_points,npeaks=1,covariance_parameterization=covariance_parameterization,return_gradient=True,return_hessian=True)

			# # print(time.time()-start)



			# dg  = peak_loss.derivative(self.fit_params)

			# start = time.time()

			# dg  = peak_loss.gradient(self.fit_params)

			# d2g = peak_loss.hessian(self.fit_params)

			# # dg  = obj.derivative(self.fit_params[nbkgr:],fit_points,npeaks=1,covariance_parameterization=covariance_parameterization,return_gradient=True,return_hessian=True)

			# # d2g = obj.hessian(self.fit_params[nbkgr:],fit_points,npeaks=1,covariance_parameterization=covariance_parameterization,return_gradient=True,return_hessian=True)

			# print(time.time()-start)



			# start = time.time()

			# # grad = numerical_gradient(self.fit_params[nbkgr:], lambda x: gaussian_mixture(x,fit_points,npeaks=1,covariance_parameterization=covariance_parameterization,return_gradient=False,return_hessian=False)[15855])

			# # hess = numerical_hessian(self.fit_params[nbkgr:],  lambda x: gaussian_mixture(x,fit_points,npeaks=1,covariance_parameterization=covariance_parameterization,return_gradient=False,return_hessian=False)[15855])

			# grad = numerical_gradient(self.fit_params, lambda x: peak_loss(x))

			# hess = numerical_hessian(self.fit_params,  lambda x: peak_loss(x))

			# print(time.time()-start)



			# np.set_printoptions(precision=1, linewidth=200)



			# print('Gradient')

			# print(dg)#[15855,:])

			# print(grad)
@@ -1316,10 +1509,10 @@ class Histogram(object): 
		data, points, edges = self.get_grid_data(return_edges=True)



		# point along each dimension

		dim_points = [points[0][:,0,0],points[1][0,:,0],points[2][0,0,:]]

		dim_points = [points[:,0,0,0],points[0,:,0,1],points[0,0,:,2]]



		# parameters of the model

		nbkgr   = 1 + self.ndims

		nbkgr   = 1 #+ self.ndims

		npeak   = self.fit_params.size - nbkgr

		ncnt    = self.ndims

		ncov    = (self.ndims*(self.ndims+1))//2
@@ -1335,7 +1528,7 @@ class Histogram(object): 
		skew     = self.fit_params[1+nbkgr+ncnt+ncov:1+nbkgr+ncnt+ncov+nskew]



		# inverse rotation matrix

		R = rotation_matrix(angles)

		R,_,_ = rotation_matrix(angles)



		# full covariance matrix

		cov_3d = R.T @ np.diag(sqrt_eig**2) @ R
@@ -1354,7 +1547,7 @@ class Histogram(object): 
			sigma_2d.append( np.sqrt(eigi) )



		# fitted model

		fit = bkgr[0]**2 + gaussian_mixture(self.fit_params[nbkgr:], points.reshape((3,-1)), covariance_parameterization='givens').reshape(data.shape)

		fit = bkgr[0]**2 + gaussian_mixture(self.fit_params[nbkgr:], points.reshape((-1,self.ndims)), covariance_parameterization='givens').reshape(data.shape)

		fit_masked = (1 if self.detector_mask is None else self.detector_mask) * fit



		# rebinned data and fit
@@ -1571,7 +1764,7 @@ class Histogram(object): 
		if bkgr_std is None: bkgr_std = peak_std + 3



		# parameters of the model

		nbkgr   = 1 + self.ndims

		nbkgr   = 1 #+ self.ndims

		npeak   = self.fit_params.size - nbkgr

		ncnt    = self.ndims

		ncov    = (self.ndims*(self.ndims+1))//2
@@ -1584,15 +1777,15 @@ class Histogram(object): 
		sqrtP    = self.fit_params[1+nbkgr+ncnt:1+nbkgr+ncnt+ncov]

		angles   = sqrtP[:nangles]

		svals    = sqrtP[nangles:]

		skew     = self.fit_params[1+nbkgr+ncnt+ncov:1+nbkgr+ncnt+ncov+nskew]

		# skew     = self.fit_params[1+nbkgr+ncnt+ncov:1+nbkgr+ncnt+ncov+nskew]





		# inverse rotation matrix

		R = rotation_matrix(angles)

		R,_,_ = rotation_matrix(angles)



		#

		data, points, edges = self.get_grid_data(return_edges=True)

		points = points.reshape((3,-1))

		points = points.reshape((-1,self.ndims))

		fit = self.fit_params[0]**2 + gaussian_mixture(self.fit_params[nbkgr:],points,npeaks=1,covariance_parameterization='givens').reshape(data.shape)



		data = data.ravel()
@@ -1671,17 +1864,19 @@ class Histogram(object): 
if __name__ == '__main__':

	run   = 43652

	loss  = 'mle'

	bins  = 32

	bins  = 64

	n_std = 4



	for peak_id in [1074]:

	# for peak_id in [1074]:

	# for peak_id in [1082]:

	# for peak_id in [1074,1082]:

	# for peak_id in [1077]:

	for peak_id in [1,6,8,11,17,22,74,1074,1077,1196,1239]:

		print(f'Processing peak {peak_id}, loss {loss}')



		hist_ws = LoadMD(f'../Result/TOPAZ_{run}_peak_{peak_id}.nxs')

		hist_ws = LoadMD(f'../Result/TOPAZ_{run}_peak_{peak_id}_{bins}.nxs')



		detector_mask = np.load(f'../Result/TOPAZ_{run}_peak_{peak_id}_mask.npy')

		detector_mask = np.load(f'../Result/TOPAZ_{run}_peak_{peak_id}_mask_{bins}.npy')

		# detector_mask = np.ones_like(hist_ws.getSignalArray())



		h = Histogram(hist_ws, detector_mask=detector_mask)
@@ -1690,7 +1885,9 @@ if __name__ == '__main__': 
		# h.fit(bins=10)



		start = time.time()

		gauss_params, fit_sucess, plot_bins = h.fit(return_bins=True)

		# gauss_params, fit_sucess, plot_bins = h.fit(return_bins=True, solver='BFGS')

		gauss_params, fit_sucess, plot_bins = h.fit_multilevel(min_level=4, return_bins=True, solver0='BFGS', solver='L-BFGS-B')

		# gauss_params, fit_sucess, plot_bins = h.fit_two_level(min_level=4, return_bins=True, solver0='Newton-CG', solver='BFGS')

		print(f"Fit: {time.time()-start} sec")



		start = time.time()
@@ -1701,6 +1898,8 @@ if __name__ == '__main__': 
		print(h.integrate())

		print(f"Integrate: {time.time()-start} sec")



		print()