rtang/PtolemyGUI
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

rewrite fit plot to add control

Browse Source
This commit is contained in:
Ryan Tang 2024-11-08 14:23:49 -05:00
parent 019877e5ea
commit def7d58217
4 changed files with 132 additions and 236 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

191
Cleopatra/FitExData.py
View File

@ -4,7 +4,7 @@ import numpy as np
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
from PyQt6.QtWidgets import ( from PyQt6.QtWidgets import (
QVBoxLayout, QWidget QVBoxLayout, QGridLayout, QWidget, QCheckBox
) )
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
@ -12,26 +12,9 @@ from matplotlib.backends.backend_qtagg import NavigationToolbar2QT as Navigation
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from ExtractXsecPy import read_DWBA from ExtractXsecPy import read_DWBA
from PlotWindow import FitPlotWindow
default_colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] #========================================================
class FitPlotWidget(QWidget):
def __init__(self, figure):
super().__init__()
self.setWindowTitle("Fit Plot")
self.resize(800, 600)
self.canvas = FigureCanvas(figure)
self.toolbar = NavigationToolbar(self.canvas, self)
layout = QVBoxLayout(self)
layout.addWidget(self.toolbar)
layout.addWidget(self.canvas)
self.setLayout(layout)
class Fitting(): class Fitting():
def __init__(self): def __init__(self):
@ -43,10 +26,16 @@ class Fitting():
self.data = [] # is a 2D array self.data = [] # is a 2D array
self.headers = [] self.headers = []
# fit parameters for a single data set
self.para = []
self.para_err = []
self.chi_squared = []
self.plot = []
def read_data(self,file_path): def read_data(self,file_path):
self.headers, self.dataX, self.data = read_DWBA(file_path) self.headers, self.dataX, self.data = read_DWBA(file_path)
self.headers = self.headers[1:]
print(self.headers)
def read_expData(self, fileName): def read_expData(self, fileName):
self.dataName_list = [] self.dataName_list = []
@ -73,8 +62,8 @@ class Fitting():
current_data = [] current_data = []
# Extract dataSet Name # Extract dataSet Name
dataName = line.split()[1] dataName = line.split()[1:]
self.dataName_list.append(dataName) self.dataName_list.append(" ".join(dataName))
# Check for fit option lines # Check for fit option lines
elif line.startswith("fit"): elif line.startswith("fit"):
@ -103,114 +92,74 @@ class Fitting():
print("Fit Options:", self.fitOption[i]) print("Fit Options:", self.fitOption[i])
print(" Data List:\n", self.expData[i]) print(" Data List:\n", self.expData[i])
def FitSingleData(self, expDataID):
print("============================================")
def FitData(self): # Get the number of fit options and cross-sections
nFit = len(self.fitOption[expDataID])
nXsec = len(self.data)
figure_list = [] # Extract experimental data (x, y, and errors)
x_exp = self.expData[expDataID][:, 0] # x positions
y_exp = self.expData[expDataID][:, 2] # y positions
y_err = self.expData[expDataID][:, 3] # y errors
for expDataID in range(len(self.expData)): self.para = []
self.para_err = []
self.chi_squared = []
for k in range(nFit):
# Get the cross-section IDs for the current fit option and strip extra spaces
xsecIDStr = self.fitOption[expDataID][k].strip()
xsecID = [int(part.strip()) for part in xsecIDStr.split('+')] if '+' in xsecIDStr else [int(xsecIDStr)]
# Ensure all cross-section IDs are valid
processFlag = True
for id in range(len(xsecID)):
if xsecID[id] >= nXsec:
print(f"Error: Requested Xsec-{xsecID[id]} exceeds the number of available cross-sections ({nXsec})")
processFlag = False
print("============================================") if processFlag == False :
continue
# Get the number of fit options and cross-sections # Define the fitting function: a weighted sum of the selected data
nFit = len(self.fitOption[expDataID]) def fit_func(x, *scale):
nXsec = len(self.data) y = np.zeros_like(x)
for p, id in enumerate(xsecID):
y += scale[p] * np.interp(x, self.dataX, self.data[id])
return y
# Extract experimental data (x, y, and errors) lower_bounds = [1e-6] * len(xsecID) # Setting a small positive lower bound
x_exp = self.expData[expDataID][:, 0] # x positions upper_bounds = [np.inf] * len(xsecID) # No upper bound
x_err = self.expData[expDataID][:, 1] # x uncertainties (errors)
y_exp = self.expData[expDataID][:, 2] # y positions
y_err = self.expData[expDataID][:, 3] # y errors
fitTheory = [] # Perform curve fitting using the fit_func and experimental data with y-errors as weights
fitTheory_lower = [] popt, pcov = curve_fit(fit_func, x_exp, y_exp, sigma=y_err, absolute_sigma=True,
fitTheory_upper = [] p0=np.ones(len(xsecID)), # Initial guess for scale parameters
bounds=(lower_bounds, upper_bounds))
para = [] self.para.append(popt)
para_err = [] perr = np.sqrt(np.diag(pcov))# Standard deviation of the parameters
chi_squared = [] self.para_err.append(perr)
for k in range(nFit): # Get the fitted model values
# Get the cross-section IDs for the current fit option and strip extra spaces y_fit = fit_func(x_exp, *popt)
xsecIDStr = self.fitOption[expDataID][k].strip() residuals = y_exp - y_fit
xsecID = [int(part.strip()) for part in xsecIDStr.split('+')] if '+' in xsecIDStr else [int(xsecIDStr)] self.chi_squared.append(np.sum((residuals / y_err) ** 2))
# Ensure all cross-section IDs are valid print(f"Fitted scale for fit {k}: {', '.join([f'{x:.3f}' for x in popt])} +/- {', '.join([f'{x:.3f}' for x in perr])} | Chi^2 : {self.chi_squared[-1]:.4f}")
processFlag = True # print(f"Fitted scale for fit {k}: {popt} +/- {perr} | Chi^2 : {chi_squared[-1]:.4f}")
for id in range(len(xsecID)):
if xsecID[id] >= nXsec:
print(f"Error: Requested Xsec-{xsecID[id]} exceeds the number of available cross-sections ({nXsec})")
processFlag = False
if processFlag == False :
continue
# Define the fitting function: a weighted sum of the selected data return self.para, self.para_err, self.chi_squared
def fit_func(x, *scale):
y = np.zeros_like(x)
for p, id in enumerate(xsecID):
y += scale[p] * np.interp(x, self.dataX, self.data[id])
return y
lower_bounds = [1e-6] * len(xsecID) # Setting a small positive lower bound
upper_bounds = [np.inf] * len(xsecID) # No upper bound
# Perform curve fitting using the fit_func and experimental data with y-errors as weights def plot_fits(self):
popt, pcov = curve_fit(fit_func, x_exp, y_exp, sigma=y_err, absolute_sigma=True,
p0=np.ones(len(xsecID)), # Initial guess for scale parameters
bounds=(lower_bounds, upper_bounds))
para.append(popt) self.plot = []
perr = np.sqrt(np.diag(pcov))# Standard deviation of the parameters
para_err.append(perr) for k , dN in enumerate(self.dataName_list):
self.FitSingleData(k)
# Get the fitted model values self.plot.append( FitPlotWindow(f"Data-{k}"))
y_fit = fit_func(x_exp, *popt) self.plot[-1].set_data(k, self.expData, self.fitOption, dN, self.dataX, self.data, self.para, self.para_err, self.chi_squared)
residuals = y_exp - y_fit self.plot[-1].plot_Fit()
chi_squared.append(np.sum((residuals / y_err) ** 2)) self.plot[-1].show()
print(f"Fitted scale for fit {k}: {', '.join([f'{x:.3f}' for x in popt])} +/- {', '.join([f'{x:.3f}' for x in perr])} | Chi^2 : {chi_squared[-1]:.4f}")
# print(f"Fitted scale for fit {k}: {popt} +/- {perr} | Chi^2 : {chi_squared[-1]:.4f}")
# Append the theoretical fit for this fit option
fitTheory.append(np.zeros_like(self.dataX))
for p, id in enumerate(xsecID):
fitTheory[-1] += popt[p] * np.interp(self.dataX, self.dataX, self.data[id])
# Optionally, you can plot the uncertainty as shaded regions (confidence intervals)
# Create the upper and lower bounds of the theoretical model with uncertainties
fitTheory_upper.append(np.zeros_like(self.dataX))
fitTheory_lower.append(np.zeros_like(self.dataX))
for p, id in enumerate(xsecID):
fitTheory_upper[-1] += (popt[p] + perr[p]) * np.interp(self.dataX, self.dataX, self.data[id])
fitTheory_lower[-1] += (popt[p] - perr[p]) * np.interp(self.dataX, self.dataX, self.data[id])
fig = plt.figure()
figure_list.append(fig)
# Plot results
plt.errorbar(x_exp, y_exp, xerr=x_err, yerr=y_err, fmt='x', label='Experimental Data', color='black', markersize = 15, elinewidth=2)
# Plot all fit theories
for i, fit in enumerate(fitTheory):
plt.plot(self.dataX, fit, label=f'Chi2:{chi_squared[i]:.3f} | Xsec:{self.fitOption[expDataID][i]}')
plt.fill_between(self.dataX, fitTheory_lower[i], fitTheory_upper[i], alpha=0.2)
# Customize plot
plt.xlabel('Angle_CM [deg]')
plt.ylabel('X-Sec [a.u.]')
plt.legend()
plt.autoscale(enable=True, axis='x', tight=True)
plt.tight_layout()
plt.yscale('log')
# Replace plt.title() with plt.text() to position the title inside the plot
plt.text(0.05, 0.05, f'Fit for Exp Data : {self.dataName_list[expDataID]}', transform=plt.gca().transAxes,
fontsize=12, verticalalignment='bottom', horizontalalignment='left', color='black')
for i, _ in enumerate(para):
plt.text(0.05, 0.1 + 0.05*i, f"Xsec-{self.fitOption[expDataID][i].strip()}: {', '.join([f'{x:.3f}' for x in para[i]])} +/- {', '.join([f'{x:.3f}' for x in para_err[i]])}" , transform=plt.gca().transAxes,
fontsize=12, verticalalignment='bottom', horizontalalignment='left', color=default_colors[i])
return figure_list

6
Cleopatra/MatPlotLibWindow.py
View File

@ -59,7 +59,7 @@ class MatPlotLibWindow(QWidget):
self.ax.plot(self.x, y, plotStyle, label=self.headers[i + 1]) self.ax.plot(self.x, y, plotStyle, label=self.headers[i + 1])
self.ax.set_xlabel("Angle_CM [Deg]") self.ax.set_xlabel("Angle_CM [Deg]")
self.ax.set_ylabel("Xsec [mb/sr]") self.ax.set_ylabel(r'$\theta_{cm}$ [deg]')
self.ax.legend(loc='upper right', frameon=True) self.ax.legend(loc='upper right', frameon=True)
# Apply log scale for y-axis if selected # Apply log scale for y-axis if selected
@ -75,7 +75,9 @@ class MatPlotLibWindow(QWidget):
self.ax.grid(False) self.ax.grid(False)
self.ax.autoscale(enable=True, axis='x', tight=True) self.ax.autoscale(enable=True, axis='x', tight=True)
self.figure.tight_layout() # self.figure.tight_layout()
self.figure.subplots_adjust(left=0.1, right=0.95, top=0.95, bottom=0.1)
self.canvas.draw_idle() self.canvas.draw_idle()

158
Cleopatra/PlotWindow.py
View File

@ -10,9 +10,10 @@ import matplotlib.pyplot as plt
from matplotlib.backends.backend_qtagg import NavigationToolbar2QT as NavigationToolbar from matplotlib.backends.backend_qtagg import NavigationToolbar2QT as NavigationToolbar
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from ExtractXsecPy import read_DWBA # Set backend to a Qt-compatible one
plt.switch_backend('QtAgg') # Or use 'Qt5Agg' if there are still issues
class PlotWindow(QWidget): class FitPlotWindow(QWidget):
def __init__(self, windowTitle): def __init__(self, windowTitle):
super().__init__() super().__init__()
@ -23,13 +24,10 @@ class PlotWindow(QWidget):
self.log_scale_checkbox = QCheckBox("Use Log Scale for Y-Axis") self.log_scale_checkbox = QCheckBox("Use Log Scale for Y-Axis")
self.log_scale_checkbox.setChecked(True) self.log_scale_checkbox.setChecked(True)
self.log_scale_checkbox.stateChanged.connect(self.plot_graph) self.log_scale_checkbox.stateChanged.connect(self.plot_Fit)
self.gridline_checkbox = QCheckBox("Show Gridlines") self.gridline_checkbox = QCheckBox("Show Gridlines")
self.gridline_checkbox.stateChanged.connect(self.plot_graph) self.gridline_checkbox.stateChanged.connect(self.plot_Fit)
self.showMarker_checkBox = QCheckBox("Show Markers")
self.showMarker_checkBox.stateChanged.connect(self.plot_graph)
self.figure, self.ax = plt.subplots() self.figure, self.ax = plt.subplots()
self.canvas = FigureCanvas(self.figure) self.canvas = FigureCanvas(self.figure)
@ -37,34 +35,13 @@ class PlotWindow(QWidget):
layout = QGridLayout(self) layout = QGridLayout(self)
layout.addWidget(self.toolbar, 0, 0, 1, 3) layout.addWidget(self.toolbar, 0, 0, 1, 3)
layout.addWidget(self.showMarker_checkBox, 1, 0) layout.addWidget(self.log_scale_checkbox, 1, 0)
layout.addWidget(self.log_scale_checkbox, 1, 1) layout.addWidget(self.gridline_checkbox, 1, 1)
layout.addWidget(self.gridline_checkbox, 1, 2)
layout.addWidget(self.canvas, 2, 0, 5, 3) layout.addWidget(self.canvas, 2, 0, 5, 3)
self.xData = [] self.setLayout(layout)
self.yData_list = []
self.header_list = []
self.yTitle = ""
self.x_exp = [] # x positions def set_data(self, ID, expData, fitOption, dataName_list, xData, yData_list, para, perr, chi_square ):
self.x_err = [] # x uncertainties (errors)
self.y_exp = [] # y positions
self.y_err = [] # y errors
self.dataName = ""
self.fitOption = []
self.para = [] # fit parameters
self.perr = [] # fit error
self.chi_square = [] # fit Chi-squared
def set_plot_data(self, xData, yData_list, header_list, yTitle):
self.xData = xData
self.yData_list = yData_list
self.header_list = header_list
self.yTitle = yTitle
def set_expData(self, expData, fitOption, dataName_list, ID):
self.x_exp = expData[ID][:, 0] self.x_exp = expData[ID][:, 0]
self.x_err = expData[ID][:, 1] self.x_err = expData[ID][:, 1]
self.y_exp = expData[ID][:, 2] self.y_exp = expData[ID][:, 2]
@ -72,37 +49,12 @@ class PlotWindow(QWidget):
self.dataName = dataName_list[ID] self.dataName = dataName_list[ID]
self.fitOption = fitOption[ID] self.fitOption = fitOption[ID]
def read_Xsec(self, file_path): self.xData = xData
headers, dataX, data = read_DWBA(file_path) self.yData_list = yData_list
self.xData = dataX
self.yData_list = data
self.header_list = headers[1:]
def set_fitResult(self, para, perr, chi_sq):
self.para = para self.para = para
self.perr = perr self.perr = perr
self.chi_square = chi_sq self.chi_square = chi_square
def plot_graph(self):
self.ax.clear()
plotStyle = '-' if not self.showMarker_checkBox.isChecked() else '-o'
for i, y in enumerate(self.yData_list):
self.ax.plot(self.xData, y, plotStyle, label=self.header_list[i])
self.ax.set_xlabel('Angle_CM [deg]')
self.ax.set_ylabel(self.yTitle)
self.ax.legend(loc='upper right', frameon=True)
# Apply log scale for y-axis if selected
if self.log_scale_checkbox.isChecked():
self.ax.set_yscale('log')
else:
self.ax.set_yscale('linear')
self.ax.autoscale(enable=True, axis='x', tight=True)
self.figure.tight_layout()
def plot_Fit(self): def plot_Fit(self):
self.ax.clear() self.ax.clear()
@ -110,55 +62,57 @@ class PlotWindow(QWidget):
self.ax.errorbar(self.x_exp, self.y_exp, xerr=self.x_err, yerr=self.y_err, self.ax.errorbar(self.x_exp, self.y_exp, xerr=self.x_err, yerr=self.y_err,
fmt='x', label='Experimental Data', color='black', markersize = 15, elinewidth=2) fmt='x', label='Experimental Data', color='black', markersize = 15, elinewidth=2)
self.ax.set_xlabel('Angle_CM [deg]') fitTheory = []
self.ax.set_ylabel(self.yTitle) fitTheory_lower = []
self.ax.legend(loc='upper right', frameon=True) fitTheory_upper = []
for k in range(len(self.fitOption)):
xsecIDStr = self.fitOption[k].strip()
xsecID = [int(part.strip()) for part in xsecIDStr.split('+')] if '+' in xsecIDStr else [int(xsecIDStr)]
fitTheory.append(np.zeros_like(self.xData))
fitTheory_upper.append(np.zeros_like(self.xData))
fitTheory_lower.append(np.zeros_like(self.xData))
for id in xsecID:
fitTheory[k] += self.para[k] * np.interp(self.xData, self.xData, self.yData_list[id])
fitTheory_upper[k] += (self.para[k] + self.perr[k]) * np.interp(self.xData, self.xData, self.yData_list[id])
fitTheory_lower[k] += (self.para[k] - self.perr[k]) * np.interp(self.xData, self.xData, self.yData_list[id])
for i, fit in enumerate(fitTheory):
self.ax.plot(self.xData, fit, label=f'Chi2:{self.chi_square[i]:.3f} | Xsec:{self.fitOption[i]}')
self.ax.fill_between(self.xData, fitTheory_lower[i], fitTheory_upper[i], alpha=0.2)
self.ax.text(0.05, 0.1 + 0.05*i, f"Xsec-{self.fitOption[i].strip()}: {', '.join([f'{x:.3f}' for x in self.para[i]])} +/- {', '.join([f'{x:.3f}' for x in self.perr[i]])}" ,
transform=plt.gca().transAxes, fontsize=12,
verticalalignment='bottom', horizontalalignment='left', color=self.default_colors[i])
# Replace plt.title() with plt.text() to position the title inside the plot
self.ax.text(0.05, 0.05, f'Fit for Exp Data : {self.dataName}', transform=plt.gca().transAxes,
fontsize=12, verticalalignment='bottom', horizontalalignment='left', color='black')
# Plot decorator
# Apply log scale for y-axis if selected # Apply log scale for y-axis if selected
if self.log_scale_checkbox.isChecked(): if self.log_scale_checkbox.isChecked():
self.ax.set_yscale('log') self.ax.set_yscale('log')
else: else:
self.ax.set_yscale('linear') self.ax.set_yscale('linear')
if self.gridline_checkbox.isChecked():
self.ax.grid(True, which='both', linestyle='--', linewidth=0.5, color='gray')
else:
self.ax.grid(False)
self.ax.set_xlabel(r'$\theta_{cm}$ [deg]')
self.ax.set_ylabel(r'd$\sigma$/d$\Omega$ [deg]')
self.ax.legend(loc='upper right', frameon=True)
self.ax.autoscale(enable=True, axis='x', tight=True) self.ax.autoscale(enable=True, axis='x', tight=True)
self.figure.tight_layout() self.figure.subplots_adjust(left=0.1, right=0.95, top=0.95, bottom=0.1)
for k in range(len(self.fitOption)):
fitTheory = []
fitTheory_lower = []
fitTheory_upper = []
xsecIDStr = self.fitOption[k].strip()
xsecID = [int(part.strip()) for part in xsecIDStr.split('+')] if '+' in xsecIDStr else [int(xsecIDStr)]
fitTheory.append(np.zeros_like(self.xData))
for p, id in enumerate(xsecID):
fitTheory += self.para[p] * np.interp(self.xData, self.xData, self.yData_list[id])
fitTheory_upper.append(np.zeros_like(self.xData))
fitTheory_lower.append(np.zeros_like(self.xData))
for p, id in enumerate(xsecID):
fitTheory_upper += (self.para[p] + self.perr[p]) * np.interp(self.xData, self.xData, self.yData_list[id])
fitTheory_lower += (self.para[p] - self.perr[p]) * np.interp(self.xData, self.xData, self.yData_list[id])
# Replace plt.title() with plt.text() to position the title inside the plot
self.ax.text(0.05, 0.05, f'Fit for Exp Data : {self.dataName}', transform=plt.gca().transAxes,
fontsize=12, verticalalignment='bottom', horizontalalignment='left', color='black')
for i, fit in enumerate(fitTheory):
self.ax.plot(self.xData, fit, label=f'Chi2:{self.chi_square[i]:.3f} | Xsec:{self.fitOption[i]}')
self.ax.fill_between(self.xData, fitTheory_lower[i], fitTheory_upper[i], alpha=0.2)
for i, _ in enumerate(self.para):
self.ax.text(0.05, 0.1 + 0.05*i, f"Xsec-{self.fitOption[i].strip()}: {', '.join([f'{x:.3f}' for x in self.para[i]])} +/- {', '.join([f'{x:.3f}' for x in self.perr[i]])}" ,
transform=plt.gca().transAxes, fontsize=12,
verticalalignment='bottom', horizontalalignment='left', color=self.default_colors[i])
self.canvas.draw_idle() self.canvas.draw_idle()

13
PyGUIQt6/PtolemyGUIPy.py
View File

@ -17,7 +17,7 @@ from CustomTextEdit import CustomTextEdit
from ExtractXsecPy import extract_xsec from ExtractXsecPy import extract_xsec
from ExWindow import ExWindow from ExWindow import ExWindow
from MatPlotLibWindow import MatPlotLibWindow from MatPlotLibWindow import MatPlotLibWindow
from FitExData import Fitting, FitPlotWidget from FitExData import Fitting
################################################## MainWindow ################################################## MainWindow
class MyWindow(QMainWindow): class MyWindow(QMainWindow):
@ -400,15 +400,9 @@ class MyWindow(QMainWindow):
def fitData(self): def fitData(self):
self.SaveExpDataFile() self.SaveExpDataFile()
self.fitCanvas = []
self.fitting.read_expData(self.ExpDataFileName) self.fitting.read_expData(self.ExpDataFileName)
self.fitting.read_data(self.DWBAFileName + ".Xsec.txt") self.fitting.read_data(self.DWBAFileName + ".Xsec.txt")
figures = self.fitting.FitData() self.fitting.plot_fits()
if figures:
for p, fig in enumerate(figures):
self.fitCanvas.append(FitPlotWidget(fig))
self.fitCanvas[-1].show()
def closeEvent(self, event): def closeEvent(self, event):
if self.plot_window: if self.plot_window:
@ -416,9 +410,6 @@ class MyWindow(QMainWindow):
if self.Ex_window: if self.Ex_window:
self.Ex_window.close() # Close the PlotWindow when MainWindow closes self.Ex_window.close() # Close the PlotWindow when MainWindow closes
self.Ex_window.__del__() self.Ex_window.__del__()
if self.fitCanvas :
for x in self.fitCanvas:
x.close()
print("============== Bye Bye ========== ") print("============== Bye Bye ========== ")
event.accept() # Accept the event to proceed with closing the main window event.accept() # Accept the event to proceed with closing the main window