From f1b2d4a580f881893ed81447fe546bda6406421d Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS_testStation" Date: Thu, 7 Nov 2024 14:05:03 -0500 Subject: [PATCH] Add FitData.py, modified python GUI for fitting. work but not perfect --- Cleopatra/CustomTextEdit.py | 27 ++++-- Cleopatra/FitExData.py | 178 ++++++++++++++++++++++++++++++++++++ PyGUIQt6/PtolemyGUIPy.py | 75 ++++++++++++++- 3 files changed, 266 insertions(+), 14 deletions(-) create mode 100644 Cleopatra/FitExData.py diff --git a/Cleopatra/CustomTextEdit.py b/Cleopatra/CustomTextEdit.py index 8925f5b..6ca29e0 100644 --- a/Cleopatra/CustomTextEdit.py +++ b/Cleopatra/CustomTextEdit.py @@ -9,17 +9,27 @@ class PythonHighlighter(QSyntaxHighlighter): super().__init__(document) # Define formatting for comments - self.comment_format = QTextCharFormat() - self.comment_format.setForeground(Qt.GlobalColor.darkGreen) + self.comment_format1 = QTextCharFormat() + self.comment_format1.setForeground(Qt.GlobalColor.darkGreen) + + self.comment_format2 = QTextCharFormat() + self.comment_format2.setForeground(Qt.GlobalColor.blue) + self.comment_format2.setBackground(Qt.GlobalColor.yellow) + + self.comment_format3 = QTextCharFormat() + self.comment_format3.setForeground(Qt.GlobalColor.magenta) def highlightBlock(self, text): - # Highlight comments - if text.startswith("#"): - self.setFormat(0, len(text), self.comment_format) + if text.startswith("#") and text.startswith("#=") == False: + self.setFormat(0, len(text), self.comment_format1) if text.startswith("$"): - self.setFormat(0, len(text), self.comment_format) + self.setFormat(0, len(text), self.comment_format1) if text.startswith("0"): - self.setFormat(0, len(text), self.comment_format) + self.setFormat(0, len(text), self.comment_format1) + if text.startswith("fit"): + self.setFormat(0, len(text), self.comment_format2) + if text.startswith("#="): + self.setFormat(0, len(text), self.comment_format3) class CustomTextEdit(QTextEdit): def __init__(self, parent=None): @@ -35,9 +45,8 @@ class CustomTextEdit(QTextEdit): # Check if Ctrl+D is pressed if event.key() == Qt.Key.Key_D and event.modifiers() == Qt.KeyboardModifier.ControlModifier: self.duplicate_line() - event.accept() # Prevent the default behavior of Ctrl+D else: - super().keyPressEvent(event) # Call the base class to handle other keys + super().keyPressEvent(event) # Call the base class to handle other keys def duplicate_line(self): cursor = self.textCursor() diff --git a/Cleopatra/FitExData.py b/Cleopatra/FitExData.py new file mode 100644 index 0000000..b6bdee3 --- /dev/null +++ b/Cleopatra/FitExData.py @@ -0,0 +1,178 @@ +#!/usr/bin/python3 + +import numpy as np +from scipy.optimize import curve_fit +import matplotlib.pyplot as plt + +from ExtractXsecPy import read_DWBA + +class Fitting: + def __init__(self): + + self.ExList = [] + self.fitOption = [] + self.expData = [] + + self.dataX = [] + self.data = [] # is a 2D array + self.headers = [] + + def read_data(self,file_path): + self.headers, self.dataX, self.data = read_DWBA(file_path) + + print(self.headers) + + def read_expData(self, fileName): + self.ExList = [] + self.fitOption = [] + self.expData = [] + + current_data = [] + + with open(fileName, "r") as file: + for line in file: + line = line.strip() + + if not line: + continue + + # Check for excitation energy lines + if line.startswith("#======================"): + # If there's an existing data block, save it + if current_data: + self.expData.append(np.array(current_data, dtype=float)) + current_data = [] + + # Extract excitation energy + Ex_energy = line.split()[1] + self.ExList.append(float(Ex_energy)) + + # Check for fit option lines + elif line.startswith("fit"): + # Add fit option parameters + fit_params = [x.strip(',') for x in line.split()[1:]] + self.fitOption.append(fit_params) + + # Parse data lines + elif not line.startswith("#"): + values = [float(x) for x in line.split()] + current_data.append(values) + + # Append the last block + if current_data: + self.expData.append(np.array(current_data, dtype=float)) + + # Convert to numpy arrays + self.ExList = np.array(self.ExList) + self.expData = [np.array(data) for data in self.expData] + + # Output the result + print("=========== Number of data set:", len(self.ExList)) + for i in range(0, len(self.ExList)): + print("-------------------------") + print(" ExList:", self.ExList[i]) + print("Fit Options:", self.fitOption[i]) + print(" Data List:\n", self.expData[i]) + + + def FitData(self): + + # Set initial offset values + x_offset, y_offset = 2000, 100 + + figure = [] + + for expDataID in range(len(self.expData)): + # Get the number of fit options and cross-sections + nFit = len(self.fitOption[expDataID]) + nXsec = len(self.data) + + # 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 + + fitTheory = [] + fitTheory_lower = [] + fitTheory_upper = [] + + 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 + + if processFlag == False : + continue + + # 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 + + + 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)) + + perr = np.sqrt(np.diag(pcov)) # Standard deviation of the parameters + print(f"Fitted scale for fit {k+1}: {popt} +/- {perr}") + + # Append the theoretical fit for this fit option + fitTheory.append(np.zeros_like(self.dataX)) + for p, id in enumerate(xsecID): + fitTheory[k] += 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[k] += (popt[p] + perr[p]) * np.interp(self.dataX, self.dataX, self.data[id]) + fitTheory_lower[k] += (popt[p] - perr[p]) * np.interp(self.dataX, self.dataX, self.data[id]) + + fig = plt.figure() + figure.append(fig) + + # Plot results + plt.errorbar(x_exp, y_exp, xerr=x_err, yerr=y_err, fmt='o', label='Experimental Data', color='blue') + + # Plot all fit theories + for i, fit in enumerate(fitTheory): + plt.plot(self.dataX, fit, label=f'Xsec:{self.fitOption[expDataID][i]} Fit') + 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.ExList[expDataID]} MeV', transform=plt.gca().transAxes, + fontsize=12, verticalalignment='bottom', horizontalalignment='left', color='black') + + manager = plt.get_current_fig_manager() + # manager.window.wm_geometry(f"+{x_offset}+{y_offset}") + manager.set_window_title(f"Exp Data : {self.ExList[expDataID]} MeV") + + x_offset += 100 + y_offset += 100 + + plt.show() \ No newline at end of file diff --git a/PyGUIQt6/PtolemyGUIPy.py b/PyGUIQt6/PtolemyGUIPy.py index e56183b..f31851c 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 ################################################## MainWindow class MyWindow(QMainWindow): @@ -25,15 +25,17 @@ class MyWindow(QMainWindow): super().__init__() self.setWindowTitle("Ptolemy GUI") - self.setGeometry(100, 100, 1000, 700) - self.setMinimumSize(400, 600) + self.setGeometry(100, 100, 1000, 800) + self.setMinimumSize(1000, 800) self.lastDWBARecord = "lastDWBA.txt" self.DWBAFileName = "" + self.ExpDataFileName = "" self.LoadLastOpenDWBASource() self.bashResult = "" self.plot_window = MatPlotLibWindow() self.Ex_window = ExWindow() + self.fitting = Fitting() # Set up Group Box for DWBA Control self.gbDWBA = QGroupBox("DWBA") @@ -146,6 +148,21 @@ class MyWindow(QMainWindow): Ex_layout.addWidget(self.sbMaXEx, 1, 1) Ex_layout.addWidget(buEx, 2, 0, 1, 2) + # ExpData Group + self.gbExpFit = QGroupBox("Exp Data Fit") + fit_layout = QVBoxLayout() + fit_layout.setAlignment(Qt.AlignmentFlag.AlignTop) + self.gbExpFit.setLayout(fit_layout) + + bnOpenExpData = QPushButton("Open ExpData") + bnOpenExpData.clicked.connect(self.LoadExpDataToTextBox) + + bnFit = QPushButton("Fit") + bnFit.clicked.connect(self.fitData) + + fit_layout.addWidget(bnOpenExpData) + fit_layout.addWidget(bnFit) + # Set up the Right Side self.bnOpenDWBASource = QPushButton("Open DWBA Source") @@ -158,6 +175,16 @@ class MyWindow(QMainWindow): self.bnSaveFile = QPushButton("Save File") self.bnSaveFile.clicked.connect(self.SaveFile) + self.bnOpenExpData = QPushButton("Open Exp Data") + self.bnOpenExpData.clicked.connect(self.OpenExpDataFile) + + self.leExpDataFileName = QLineEdit("") + self.leExpDataFileName.setReadOnly(True) + self.leExpDataFileName.setText(self.ExpDataFileName) + + self.bnSaveExpDataFile = QPushButton("Save Exp Data File") + self.bnSaveExpDataFile.clicked.connect(self.SaveExpDataFile) + self.text_edit = CustomTextEdit(self) self.leStatus = QLineEdit("") @@ -165,17 +192,25 @@ class MyWindow(QMainWindow): self.LoadFileToTextBox(self.DWBAFileName, True) + self.bnSaveExpDataFile.setEnabled(False) + # Set up the layout layout = QGridLayout() # layout.addWidget(self.gbDWBA, 0, 0, 7, 1) layout.addWidget(self.gbDWBA, 0, 0, 5, 1) layout.addWidget(self.gbEx, 5, 0, 2, 1) + layout.addWidget(self.gbExpFit, 7, 0, 2, 1) layout.addWidget(self.bnOpenDWBASource, 0, 1) layout.addWidget(self.leFileName, 0, 2, 1, 3) layout.addWidget(self.bnSaveFile, 0, 5) - layout.addWidget(self.text_edit, 1, 1, 5, 5) - layout.addWidget(self.leStatus, 6, 1, 1, 5) + + layout.addWidget(self.bnOpenExpData, 1, 1) + layout.addWidget(self.leExpDataFileName, 1, 2, 1, 3) + layout.addWidget(self.bnSaveExpDataFile, 1, 5) + + layout.addWidget(self.text_edit, 2, 1, 6, 5) + layout.addWidget(self.leStatus, 8, 1, 1, 5) for i in range(layout.columnCount()) : layout.setColumnStretch(i, 1) @@ -192,12 +227,16 @@ class MyWindow(QMainWindow): try : with open(self.lastDWBARecord, 'r') as file: self.DWBAFileName = file.readline().strip() + self.ExpDataFileName = file.readline().strip() except: self.DWBAFileName = "DWBA" + self.ExpDataFileName = "" def SaveLastOpenDWBASource(self): with open(self.lastDWBARecord, 'w') as file: file.write(self.DWBAFileName) + file.write("\n") + file.write(self.ExpDataFileName) def OnOffXsecOption(self): self.cbXsec.setEnabled(self.chkExtracrXsec.isChecked()) @@ -209,6 +248,20 @@ class MyWindow(QMainWindow): self.leFileName.setText(self.DWBAFileName) self.LoadFileToTextBox(self.DWBAFileName) self.SaveLastOpenDWBASource() + self.bnSaveExpDataFile.setEnabled(False) + + def OpenExpDataFile(self): + file_path, _ = QFileDialog.getOpenFileName(self, "Open File", "", "Text File (*.txt)") + if file_path: + self.ExpDataFileName = file_path + self.leExpDataFileName.setText(self.ExpDataFileName) + self.LoadFileToTextBox(self.ExpDataFileName) + self.bnSaveExpDataFile.setEnabled(True) + self.SaveLastOpenDWBASource() + + def LoadExpDataToTextBox(self): + self.LoadFileToTextBox(self.ExpDataFileName) + self.bnSaveExpDataFile.setEnabled(True) def LoadFileToTextBox(self, fileName, moveToButton = False): # print(fileName) @@ -227,12 +280,20 @@ class MyWindow(QMainWindow): except Exception as e: self.text_edit.setText(f"Failed to load file:\n{e}") self.leStatus.setText(f"Failed to load file:\n{e}") + + self.bnSaveExpDataFile.setEnabled(False) def SaveFile(self): file_path = self.leFileName.text() with open(file_path, 'w') as file: file.write(self.text_edit.toPlainText()) self.leStatus.setText(f"File saved to: {file_path}") + + def SaveExpDataFile(self): + file_path = self.leExpDataFileName.text() + with open(file_path, 'w') as file: + file.write(self.text_edit.toPlainText()) + self.leStatus.setText(f"File saved to: {file_path}") def DeleteinOutXsecFiles(self): if os.path.exists(self.DWBAFileName + ".in"): @@ -303,6 +364,10 @@ class MyWindow(QMainWindow): self.Ex_window.plot_Ex_graph() self.Ex_window.show() + def fitData(self): + self.fitting.read_expData(self.ExpDataFileName) + self.fitting.read_data(self.DWBAFileName + ".Xsec.txt") + self.fitting.FitData() def closeEvent(self, event): if self.plot_window: