diff --git a/Cleopatra/FitExData.py b/Cleopatra/FitExData.py index 9f5a469..614cf9c 100644 --- a/Cleopatra/FitExData.py +++ b/Cleopatra/FitExData.py @@ -4,7 +4,7 @@ import numpy as np from scipy.optimize import curve_fit from PyQt6.QtWidgets import ( - QVBoxLayout, QWidget + QVBoxLayout, QGridLayout, QWidget, QCheckBox ) 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 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(): def __init__(self): @@ -43,10 +26,16 @@ class Fitting(): self.data = [] # is a 2D array 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): self.headers, self.dataX, self.data = read_DWBA(file_path) - - print(self.headers) + self.headers = self.headers[1:] def read_expData(self, fileName): self.dataName_list = [] @@ -73,8 +62,8 @@ class Fitting(): current_data = [] # Extract dataSet Name - dataName = line.split()[1] - self.dataName_list.append(dataName) + dataName = line.split()[1:] + self.dataName_list.append(" ".join(dataName)) # Check for fit option lines elif line.startswith("fit"): @@ -103,114 +92,74 @@ class Fitting(): print("Fit Options:", self.fitOption[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 - nFit = len(self.fitOption[expDataID]) - nXsec = len(self.data) + # Define the fitting function: a weighted sum of the selected data + 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 - # Extract experimental data (x, y, and errors) - x_exp = self.expData[expDataID][:, 0] # x positions - 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 + lower_bounds = [1e-6] * len(xsecID) # Setting a small positive lower bound + upper_bounds = [np.inf] * len(xsecID) # No upper bound - fitTheory = [] - fitTheory_lower = [] - fitTheory_upper = [] + # Perform curve fitting using the fit_func and experimental data with y-errors as weights + 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 = [] - para_err = [] - chi_squared = [] + self.para.append(popt) + perr = np.sqrt(np.diag(pcov))# Standard deviation of the parameters + self.para_err.append(perr) - 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)] + # Get the fitted model values + y_fit = fit_func(x_exp, *popt) + residuals = y_exp - y_fit + self.chi_squared.append(np.sum((residuals / y_err) ** 2)) - # 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 - - if processFlag == False : - continue + 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}") + # print(f"Fitted scale for fit {k}: {popt} +/- {perr} | Chi^2 : {chi_squared[-1]:.4f}") - # Define the fitting function: a weighted sum of the selected data - 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 + return self.para, self.para_err, self.chi_squared - 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 - 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)) + def plot_fits(self): - para.append(popt) - perr = np.sqrt(np.diag(pcov))# Standard deviation of the parameters - para_err.append(perr) - - # Get the fitted model values - y_fit = fit_func(x_exp, *popt) - residuals = y_exp - y_fit - chi_squared.append(np.sum((residuals / y_err) ** 2)) - - 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 \ No newline at end of file + self.plot = [] + + for k , dN in enumerate(self.dataName_list): + self.FitSingleData(k) + self.plot.append( FitPlotWindow(f"Data-{k}")) + self.plot[-1].set_data(k, self.expData, self.fitOption, dN, self.dataX, self.data, self.para, self.para_err, self.chi_squared) + self.plot[-1].plot_Fit() + self.plot[-1].show() diff --git a/Cleopatra/MatPlotLibWindow.py b/Cleopatra/MatPlotLibWindow.py index 8eb1376..b77bee2 100644 --- a/Cleopatra/MatPlotLibWindow.py +++ b/Cleopatra/MatPlotLibWindow.py @@ -59,7 +59,7 @@ class MatPlotLibWindow(QWidget): self.ax.plot(self.x, y, plotStyle, label=self.headers[i + 1]) 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) # Apply log scale for y-axis if selected @@ -75,7 +75,9 @@ class MatPlotLibWindow(QWidget): self.ax.grid(False) 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() diff --git a/Cleopatra/PlotWindow.py b/Cleopatra/PlotWindow.py index e54bf5e..cd6b936 100644 --- a/Cleopatra/PlotWindow.py +++ b/Cleopatra/PlotWindow.py @@ -10,9 +10,10 @@ import matplotlib.pyplot as plt from matplotlib.backends.backend_qtagg import NavigationToolbar2QT as NavigationToolbar 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): super().__init__() @@ -23,13 +24,10 @@ class PlotWindow(QWidget): self.log_scale_checkbox = QCheckBox("Use Log Scale for Y-Axis") 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.stateChanged.connect(self.plot_graph) - - self.showMarker_checkBox = QCheckBox("Show Markers") - self.showMarker_checkBox.stateChanged.connect(self.plot_graph) + self.gridline_checkbox.stateChanged.connect(self.plot_Fit) self.figure, self.ax = plt.subplots() self.canvas = FigureCanvas(self.figure) @@ -37,34 +35,13 @@ class PlotWindow(QWidget): layout = QGridLayout(self) layout.addWidget(self.toolbar, 0, 0, 1, 3) - layout.addWidget(self.showMarker_checkBox, 1, 0) - layout.addWidget(self.log_scale_checkbox, 1, 1) - layout.addWidget(self.gridline_checkbox, 1, 2) + layout.addWidget(self.log_scale_checkbox, 1, 0) + layout.addWidget(self.gridline_checkbox, 1, 1) layout.addWidget(self.canvas, 2, 0, 5, 3) - self.xData = [] - self.yData_list = [] - self.header_list = [] - self.yTitle = "" + self.setLayout(layout) - self.x_exp = [] # x positions - 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): + def set_data(self, ID, expData, fitOption, dataName_list, xData, yData_list, para, perr, chi_square ): self.x_exp = expData[ID][:, 0] self.x_err = expData[ID][:, 1] self.y_exp = expData[ID][:, 2] @@ -72,37 +49,12 @@ class PlotWindow(QWidget): self.dataName = dataName_list[ID] self.fitOption = fitOption[ID] - def read_Xsec(self, file_path): - headers, dataX, data = read_DWBA(file_path) - self.xData = dataX - self.yData_list = data - self.header_list = headers[1:] + self.xData = xData + self.yData_list = yData_list - def set_fitResult(self, para, perr, chi_sq): self.para = para self.perr = perr - self.chi_square = chi_sq - - 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() + self.chi_square = chi_square def plot_Fit(self): 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, fmt='x', label='Experimental Data', color='black', markersize = 15, elinewidth=2) - self.ax.set_xlabel('Angle_CM [deg]') - self.ax.set_ylabel(self.yTitle) - self.ax.legend(loc='upper right', frameon=True) + fitTheory = [] + fitTheory_lower = [] + 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 if self.log_scale_checkbox.isChecked(): self.ax.set_yscale('log') else: 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.figure.tight_layout() - - 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.figure.subplots_adjust(left=0.1, right=0.95, top=0.95, bottom=0.1) self.canvas.draw_idle() - - - - - - diff --git a/PyGUIQt6/PtolemyGUIPy.py b/PyGUIQt6/PtolemyGUIPy.py index 506bef7..10019b0 100755 --- a/PyGUIQt6/PtolemyGUIPy.py +++ b/PyGUIQt6/PtolemyGUIPy.py @@ -17,7 +17,7 @@ from CustomTextEdit import CustomTextEdit from ExtractXsecPy import extract_xsec from ExWindow import ExWindow from MatPlotLibWindow import MatPlotLibWindow -from FitExData import Fitting, FitPlotWidget +from FitExData import Fitting ################################################## MainWindow class MyWindow(QMainWindow): @@ -400,15 +400,9 @@ class MyWindow(QMainWindow): def fitData(self): self.SaveExpDataFile() - self.fitCanvas = [] self.fitting.read_expData(self.ExpDataFileName) self.fitting.read_data(self.DWBAFileName + ".Xsec.txt") - figures = self.fitting.FitData() - - if figures: - for p, fig in enumerate(figures): - self.fitCanvas.append(FitPlotWidget(fig)) - self.fitCanvas[-1].show() + self.fitting.plot_fits() def closeEvent(self, event): if self.plot_window: @@ -416,9 +410,6 @@ class MyWindow(QMainWindow): if self.Ex_window: self.Ex_window.close() # Close the PlotWindow when MainWindow closes self.Ex_window.__del__() - if self.fitCanvas : - for x in self.fitCanvas: - x.close() print("============== Bye Bye ========== ") event.accept() # Accept the event to proceed with closing the main window