Commit 1db49332 authored by Viktor Reshniak's avatar Viktor Reshniak
Browse files

fits

parent e564e9c7

peakintegration/peakfitter.py

+91 −91
Original line number Diff line number Diff line
@@ -1548,16 +1548,16 @@ class PeaksIntegrator(object): 
					params_ubnd = [ np.inf for _ in range(2)] + [c+peak_rad/5 for c in ellipsoid[:3]] + [ np.inf for phi in ellipsoid[3:6]] + [min(2.0*ax,0.9*dst2nn/peak_std,0.5*box_size/peak_std) for ax in ellipsoid[6:]]

					params_reg  = None



					# box_size has been changed, need to rebin anyway

					binned_peak = None

				else:

					# ellipsoid parameters with confidence and prediction intervals

					ellparams, params_confint, params_predint = ellipsoid_predictor(np.array([self.peaks[peak_id].getWavelength()]))

					ellparams, params_confint, params_predint = ellparams[0], params_confint[0], params_predint[0]

					ellparams, params_confint, params_predint = ellipsoid_predictor(self.peak_wavelen[peak_id])



					# update box size

					peak_rad = np.max(ellparams[-3:])*peak_std

					if peak_rad>5*np.min(ellparams[-3:])*peak_std:

						# adjust for too 'flat' peaks

						# adjust for too 'flat' peaks - replace box_std with peak_std+2

						box_size = self.valid_box_size(peak_id, 2*((peak_std+2)/peak_std)*peak_rad)

					else:

						box_size = self.valid_box_size(peak_id, 2*(box_std/peak_std)*peak_rad)
@@ -1574,17 +1574,15 @@ class PeaksIntegrator(object): 


					ellipsoid = [c for c in self.peak_centers[peak_id]] + list(ellparams)





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

					# parameter bounds

					# predict intensity



					# predicted intensity

					binned_peak = self.bin_peak(peak_id, bins, hkl_box_size=box_size, smooth=False)

					points,_    = histogram_grid(binned_peak)

					if self.peak_funs[peak_id] is None:

						self.peak_funs[peak_id] = Gaussian(3, 'givens')

					self.peak_funs[peak_id].set_ellipsoid(np.array(ellipsoid[:3]), np.array(ellipsoid[3:]))

					peak_mask = self.peak_funs[peak_id].is_inside_ellipsoid(std=_peak_std, x=points.reshape((-1,3))).reshape([bins]*3)

					I,sig = integrate_peak(binned_peak.getSignalArray().copy(), peak_mask, ~peak_mask, valid_mask=None)

					peak_mask   = Gaussian(3,'givens').is_inside_ellipsoid(std=_peak_std, x=points.reshape((-1,3)), params=[1]+ellipsoid).reshape([bins]*3)



					I,sig = integrate_peak(binned_peak.getSignalArray(), peak_mask, ~peak_mask, valid_mask=None)

					sigI  = sig / max(I,0.1)



					# do not fit weak peaks - use predictor
@@ -1592,7 +1590,9 @@ class PeaksIntegrator(object): 
						self.peaks[peak_id].setIntensity(I)

						self.peaks[peak_id].setSigmaIntensity(sig)

						return

						# continue



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

					# parameter bounds



					resolution = box_size/bins
@@ -1606,16 +1606,15 @@ class PeaksIntegrator(object): 
					# params_reg  = [ 0      for _ in range(2)] + [0             for c in ellipsoid[:3]] + [1/p for p in params_confint[:3]] + [1]*3



			else:

				# id of main peak

				main_id = self.main_ids[peak_id]



				if ellipsoid_predictor is None:

					# id of the valid main peak

				if main_id is None or self.peak_funs[main_id] is None:

						for main_id in self._peaks_nn_ind['centers'][peak_id,:]:

					for main_id in self._peaks_nn_ind['wavelengths'][peak_id,:]:

						if not self.is_satellite(main_id) and self.peak_funs[main_id] is not None:

							break



					# estimate box size by looking at satellites of the main peak, choose the largest

				if ellipsoid_predictor is None:

					# # estimate box size by looking at satellites of the main peak, choose the largest

					# dst2main = max([self.dst_between_peaks(pid, main_id) for pid in self.sat_ids[main_id]])

					box_size = self.valid_box_size(peak_id, 2*0.8*dst2main)
@@ -1623,51 +1622,55 @@ class PeaksIntegrator(object): 
					main_rad = max(self.peak_funs[main_id].semiaxes) * peak_std

					max_rad  = dst2main - main_rad



					# histogram scale averaging

					# smooth_op = self.smooth_op(self.sat_ids[main_id],box_size,bins)

					smooth_op = self.smooth_op(peak_id,box_size,bins)



					datas = []

					# for pid in self.sat_ids[main_id]:

					for pid in [peak_id]:

						points,_  = histogram_grid(self.bin_peak(pid, bins, box_size, smooth=False))

						main_mask = self.peak_funs[main_id].is_inside_ellipsoid(std=_peak_std, x=points.reshape((-1,3))).reshape([bins]*3)

						binned_peak   = self.bin_peak(pid, bins, box_size, smooth=False)

						binned_points = histogram_grid(binned_peak)[0]

						binned_data   = binned_peak.getSignalArray().copy()



						# mask out main peak

						main_mask = self.peak_funs[main_id].is_inside_ellipsoid(std=_peak_std, x=binned_points.reshape((-1,3))).reshape([bins]*3)

						binned_data[main_mask] = 0

						datas.append(binned_data)

					binned_data = sum(datas)



						data = self.bin_peak(pid, bins, hkl_box_size=box_size, smooth=False).getSignalArray().copy()

						data[main_mask] = 0

						datas.append(data)

					data = sum(datas)

					# # histogram scale averaging

					# smooth_op = self.smooth_op(self.sat_ids[main_id],box_size,bins)

					smooth_op   = self.smooth_op(peak_id,box_size,bins,binned_hist=binned_data)

					smooth_data = smooth_op(binned_data.ravel()).reshape(binned_data.shape)



					if smooth_op(data.ravel()).max()<10:

					if smooth_data.max()<weak_density_thresh:

						return

						# continue



					# plot_volumes([smooth_op(d.ravel()).reshape([bins]*3) for d in datas]+[smooth_op(data.ravel()).reshape([bins]*3)], plot_axes=False, show=True, elevation=15, azimuth=-40)



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



					ellipsoid = isolate_peak(points=points, data=smooth_op(data.ravel()).reshape(data.shape), nthresh=20, iters=2)

					ellipsoid = isolate_peak(points=binned_points, data=smooth_data, nthresh=20, iters=2)

					# ellipsoid can't exceed max_rad

					for i in range(3,6):

						ellipsoid[i] = min(ellipsoid[i], 0.9*max_rad)

					peak_rad = np.max(ellipsoid[-3:])

					if peak_rad>5*np.min(ellipsoid[-3:]):

						# adjust for too 'flat' peaks

						# adjust for too 'flat' peaks - replace box_std with peak_std+2

						box_size = self.valid_box_size(peak_id, 2*((peak_std+2)/peak_std)*peak_rad)

					else:

						box_size = self.valid_box_size(peak_id, 2*(box_std/peak_std)*peak_rad)



					# scale semiaxes to 1 std

					# rescale semiaxes to 1 std

					ellipsoid[-3:] /= peak_std



					# center at box center

					ellipsoid = [c for c in self.peak_centers[peak_id]] + list(ellipsoid[3:])



					# reset angles to main peak angles

					ellipsoid[3:6] = self.peak_funs[main_id].angles



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



					# 'fake' data of combined satellites

					binned_peak = self.bin_peak(peak_id, bins, box_size, smooth=True)

					# binned_peak = self.bin_peak(peak_id, bins, box_size, smooth=True)

					# binned_peak.setSignalArray(data)

					# plot_volumes([binned_peak.getSignalArray()], plot_axes=False, show=True, elevation=15, azimuth=-40)
@@ -1680,10 +1683,11 @@ class PeaksIntegrator(object): 
					params_ubnd = [ np.inf for _ in range(2)] + [c+max_rad/5 for c in ellipsoid[:3]] + [phi+np.pi/10 for phi in ellipsoid[3:6]] + [min(1.3*ax,0.9*max_rad/peak_std,box_size/2/peak_std) for ax in ellipsoid[6:]]

					params_reg  = None



					# box_size has been changed, need to rebin anyway

					binned_peak = None

				else:

					# ellipsoid parameters of main peak with confidence and prediction intervals

					ellparams, params_confint, params_predint = ellipsoid_predictor(np.array([self.peaks[peak_id].getWavelength()]))

					ellparams, params_confint, params_predint = ellparams[0], params_confint[0], params_predint[0]

					ellparams, params_confint, params_predint = ellipsoid_predictor(self.peak_wavelen[peak_id])



					# # estimate box size by looking at peak neighbors

					# box_size = self.valid_box_size(peak_id, 2*0.8*dst2main)
@@ -1696,28 +1700,23 @@ class PeaksIntegrator(object): 
					else:

						box_size = self.valid_box_size(peak_id, 2*(box_std/peak_std)*peak_rad)



					# # isolate peak and initialize ellipsoid

					# binned_peak = self.bin_peak(peak_id, 64, hkl_box_size=box_size, smooth=False)

					# if binned_peak.getNEvents()<=50:

					# 	continue



					# main_binned = self.bin_peak(self.main_ids[peak_id], 64, hkl_box_size=0.4, smooth=False)

					# main_data = main_binned.getSignalArray().copy()



					ellipsoid = [c for c in self.peak_centers[peak_id]] + list(ellparams)



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

					# predict intensity



					binned_peak = self.bin_peak(peak_id, bins, hkl_box_size=box_size, smooth=False)

					points,_    = histogram_grid(binned_peak)

					if self.peak_funs[peak_id] is None:

						self.peak_funs[peak_id] = Gaussian(3, 'givens')

					self.peak_funs[peak_id].set_ellipsoid(np.array(ellipsoid[:3]), np.array(ellipsoid[3:]))

					peak_mask = self.peak_funs[peak_id].is_inside_ellipsoid(std=_peak_std, x=points.reshape((-1,3))).reshape([bins]*3)

					bkgr_mask = ~peak_mask

					peak_mask   = Gaussian(3,'givens').is_inside_ellipsoid(std=_peak_std, x=points.reshape((-1,3)), params=[1]+ellipsoid).reshape([bins]*3)

					if main_id is not None:

						main_mask = self.peak_funs[main_id].is_inside_ellipsoid(std=_peak_std, x=points.reshape((-1,3))).reshape([bins]*3)

					else:

						main_mask = np.full_like(peak_mask, fill_value=False, dtype=bool)

					I, sig = integrate_peak(binned_peak.getSignalArray().copy(), peak_mask, bkgr_mask, valid_mask=~main_mask)



					I,sig = integrate_peak(binned_peak.getSignalArray(), peak_mask, ~peak_mask, valid_mask=~main_mask)

					sigI  = sig / max(I,0.1)



					# do not fit weak peaks - use predictor
@@ -1725,17 +1724,42 @@ class PeaksIntegrator(object): 
						self.peaks[peak_id].setIntensity(I)

						self.peaks[peak_id].setSigmaIntensity(sig)

						return

						# continue



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



					resolution = box_size/bins

					# binned_peak = self.bin_peak(peak_id, bins, hkl_box_size=box_size, smooth=False)

					# points,_   = histogram_grid(binned_peak)

					# if self.peak_funs[peak_id] is None:

					# 	self.peak_funs[peak_id] = Gaussian(3, 'givens')

					# self.peak_funs[peak_id].set_ellipsoid(np.array(ellipsoid[:3]), np.array(ellipsoid[3:]))

					# peak_mask = self.peak_funs[peak_id].is_inside_ellipsoid(std=_peak_std, x=points.reshape((-1,3))).reshape([bins]*3)

					# bkgr_mask = ~peak_mask

					# if main_id is not None:

					# 	main_mask = self.peak_funs[main_id].is_inside_ellipsoid(std=_peak_std, x=points.reshape((-1,3))).reshape([bins]*3)

					# else:

					# 	main_mask = np.full_like(peak_mask, fill_value=False, dtype=bool)

					# I, sig = integrate_peak(binned_peak.getSignalArray().copy(), peak_mask, bkgr_mask, valid_mask=~main_mask)

					# sigI = sig / max(I,0.1)



					# plot_volumes([binned_peak.getSignalArray(),self.smooth_op(peak_id,box_size,bins)(binned_peak.getSignalArray().ravel()).reshape([bins]*3)], plot_axes=False, show=True, elevation=15, azimuth=-40)



					# # do not fit weak peaks - use predictor

					# if I<50 and len(peak_ids)>1:

					# 	self.peaks[peak_id].setIntensity(I)

					# 	self.peaks[peak_id].setSigmaIntensity(sig)

					# 	return

					# 	# continue



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

					# parameter bounds



					resolution = box_size/bins



					#                   intensities                                center                                                  angle                                         axes

					params_init = [ None   for _ in range(2)] + [c            for c in ellipsoid[:3]] + [phi          for phi in ellipsoid[3:6]] + [ax                                                    for ax in ellipsoid[6:]]

					params_lbnd = [-np.inf for _ in range(2)] + [c-peak_rad/5 for c in ellipsoid[:3]] + [phi-np.pi/20 for phi in ellipsoid[3:6]] + [max(0.2*ax,resolution/peak_std)                       for ax in ellipsoid[6:]]

					params_ubnd = [ np.inf for _ in range(2)] + [c+peak_rad/5 for c in ellipsoid[:3]] + [phi+np.pi/20 for phi in ellipsoid[3:6]] + [min(1.3*ax,0.9*dst2nn/peak_std,box_size/2.5/peak_std) for ax in ellipsoid[6:]]

					params_reg  = [ 0      for _ in range(2)] + [0            for c in ellipsoid[:3]] + [sigI**2/p**2 for p in params_predint]

					params_reg  = [ 0      for _ in range(2)] + [0            for c in ellipsoid[:3]] + [min(sigI**2/p**2,1.e5) for p in params_predint]

					# params_lbnd = [-np.inf for _ in range(2)] + [c-box_size/10 for c in ellipsoid[:3]] + [-np.inf for phi in ellipsoid[3:6]] + [max(0.2*ax,box_size/10/peak_std)  for ax in ellipsoid[6:]]

					# params_ubnd = [ np.inf for _ in range(2)] + [c+box_size/10 for c in ellipsoid[:3]] + [ np.inf for phi in ellipsoid[3:6]] + [min(1.3*ax,box_size/2.5/peak_std) for ax in ellipsoid[6:]]

					# params_lbnd = [-np.inf for _ in range(2)] + [c-box_size/10 for c in self.peak_centers[peak_id]] + [phi-np.pi/10 for phi in ellparams[:3]] + [max(ax-box_size/10,box_size/6/peak_std)   for ax in ellipsoid[-3:]]
@@ -1749,37 +1773,13 @@ class PeaksIntegrator(object): 
			# attempt fitting the peak



			h = self.fit_peak(peak_id, bins, box_size, binned_peak=binned_peak, params_init=params_init, params_lbnd=params_lbnd, params_ubnd=params_ubnd, ellipsoid=ellipsoid, peak_rad=peak_rad, solver='BFGS', tol=tol, multilevel=multilevel, params_reg=params_reg)

			# if ellipsoid_predictor is not None:

			# 	h = self.fit_peak(peak_id, bins, box_size, params_init=params_init, params_lbnd=params_lbnd, params_ubnd=params_ubnd, ellipsoid=ellipsoid, peak_rad=peak_rad, solver='BFGS', tol=tol, multilevel=multilevel, params_reg=params_reg)

			# else:

			# 	params_lbnd = [-np.inf for _ in range(2)] + params_init[2:]

			# 	params_ubnd = [ np.inf for _ in range(2)] + params_init[2:]

			# 	h = self.fit_peak(peak_id, bins, box_size, params_init=params_init, params_lbnd=params_lbnd, params_ubnd=params_ubnd, ellipsoid=ellipsoid, peak_rad=peak_rad, solver='BFGS', tol=tol, multilevel=multilevel, params_reg=params_reg)



			# if ellipsoid_predictor is not None and (h.status_code!=0):

			# 	ellipsoid = [c for c in self.peak_centers[peak_id]] + [0]*3 + [box_size/4/peak_std]*3

			# 	if self.peak_mnps[peak_id]=='0,0,0':

			# 		params_init = [ None for _ in range(2)] + [c for c in self.peak_centers[peak_id]] + [phi for phi in ellparams[:3]] + [ax for ax in ellparams[3:]]

			# 	else:

			# 		params_init = [ None for _ in range(2)] + ellipsoid

			# 	params_lbnd = [-np.inf for _ in range(2)] + [c-box_size/20 for c in self.peak_centers[peak_id]] + params_init[5:]

			# 	params_ubnd = [ np.inf for _ in range(2)] + [c+box_size/20 for c in self.peak_centers[peak_id]] + params_init[5:]

			# 	params_reg  = [0]*2 + [1]*3 + [1.e6]*6

			# 	h = self.fit_peak(peak_id, bins, box_size, params_init=params_init, params_lbnd=params_lbnd, params_ubnd=params_ubnd, ellipsoid=ellipsoid, peak_rad=peak_rad, solver='BFGS', tol=tol, multilevel=multilevel, params_reg=params_reg)





			# if h.status_code==1:

			# 	h = self.fit_peak(peak_id, bins, box_size, solver='L-BFGS-B')

			# elif h.status_code==2:

			# 	h = self.fit_peak(peak_id, bins, box_size, solver='L-BFGS-B')#, params_init=h.fit_params)





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



			self.valid_masks[peak_id] = h.fit_mask



			self.integrate_peak(peak_id, h)





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



			# if _verbose:
@@ -1792,17 +1792,13 @@ class PeaksIntegrator(object): 
			if plot_path is not None:

				prefix = f'run{self.peaks_ws.getRunNumber()}_{self.peak_banks[peak_id]}_peak{peak_id}_wavel{self.peak_wavelen[peak_id]:.1f}_bins{bins}_box{box_size:.2f}_I{int(self.peaks[peak_id].getIntensity())}_status{h.status_code}'

				if len(peak_ids)==1:

					h.plot(peak_std=peak_std, plot_path=Path(plot_path).parent, bkgr_std=10, prefix=prefix, plot_mask='fit', show_empty=True, evolution=False)

					h.plot(peak_std=peak_std, plot_path=Path(plot_path).parent, bkgr_std=_bkgr_std, prefix=prefix, plot_mask='fit', show_empty=True, evolution=False)

				else:

					h.plot(peak_std=peak_std, plot_path=Path(plot_path,f'mnp_{self.peak_mnps[peak_id]}',f'{h.status_message}'), bkgr_std=10, prefix=prefix, evolution=False)

				# elif ellipsoid_predictor is not None:

				# 	h.plot(peak_std=peak_std, plot_path=Path(plot_path,'predicted',f'{h.status_message}'), bkgr_std=10, prefix=prefix, evolution=False)

					h.plot(peak_std=peak_std, plot_path=Path(plot_path,f'mnp_{self.peak_mnps[peak_id]}',f'{h.status_message}'), bkgr_std=_bkgr_std, prefix=prefix, evolution=False)



			# self.plot_peak(peak_id, bins=[16,32], hkl_box_size=box_size, plot_path=f'data', prefix=f'peak_{peak_id}_boxsize_{box_size:.3f}_{prev_box_size:.3f}_status_{h.status_code}')









		for peak_id in peak_ids:

			try_catch(peak_id)

			# try:
@@ -1843,7 +1839,7 @@ class PeaksIntegrator(object): 
		# determine order of peaks

		if sort_peaks:

			# counts of peak boxes

			box_size = min(self.nearest_hkl_neighbors[0].min(), self.max_hkl_box_size) / 3

			box_size = min(self.nearest_hkl_dst.min(), self.max_hkl_box_size) / 3

			counts = np.zeros((self.npeaks,), dtype=int)

			for peak_id in range(self.npeaks):

				counts[peak_id] = self.bin_peak(peak_id,bins=1,hkl_box_size=box_size).getNEvents()
@@ -1887,6 +1883,7 @@ class PeaksIntegrator(object): 
		# 	ellipsoid_predictor = None

		# exit()



		if len(peak_ids_sat)>0:

			ellipsoid_predictor_sat = self.fit_ellipsoid(peak_ids_sat, intensity_thresh=100, plot_path=make_path(plot_path,f'iter1_sat_ellipsoid_predictor.png'), base_predictor=ellipsoid_predictor_main)

		# exit()
@@ -1895,12 +1892,15 @@ class PeaksIntegrator(object): 
		# second iteration



		self.fit_peak_ids(peak_ids_main, bins, ellipsoid_predictor=ellipsoid_predictor_main, plot_path=make_path(plot_path,f'iter_{2}'))

		if len(peak_ids_sat)>0:

			self.fit_peak_ids(peak_ids_sat, bins, ellipsoid_predictor=ellipsoid_predictor_sat, plot_path=make_path(plot_path,f'iter_{2}'))

		self.save_parameters_of_integrated_peaks(peak_ids, save_path=make_path(self.save_path,f'params_iter_{2}'))

		# self.load_parameters_of_integrated_peaks(load_path=make_path(self.save_path,f'params_iter_{2}'))



		self.integrate_loaded(bins)



		ellipsoid_predictor_main = self.fit_ellipsoid(peak_ids_main, intensity_thresh=strong_counts_thresh, plot_path=make_path(plot_path,f'iter2_main_ellipsoid_predictor.png'))

		if len(peak_ids_sat)>0:

			ellipsoid_predictor_sat  = self.fit_ellipsoid(peak_ids_sat,  intensity_thresh=10,                   plot_path=make_path(plot_path,f'iter2_sat_ellipsoid_predictor.png'))

		# exit()