Create log file vesuvio

src/mvesuvio/analysis_fitting.py

@@ -9,8 +9,11 @@
 import jacobi
 import time
 import re
+from mantid.kernel import logger
+from mantid.simpleapi import AnalysisDataService
 
 from mvesuvio.util import handle_config
+from mvesuvio.util.analysis_helpers import print_table_workspace, pass_data_into_ws
 
 repoPath = Path(__file__).absolute().parent  # Path to the repository
 
@@ -20,7 +23,7 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF):
     try:
         wsTOF = mtd[wsTOF]
     except KeyError: 
-        print(f"Workspace to fit {wsTOF} not found.")
+        logger.notice(f"Workspace to fit {wsTOF} not found.")
         return
 
     wsResSum, wsRes = calculateMantidResolutionFirstMass(IC, yFitIC, wsTOF)
@@ -39,7 +42,7 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF):
     fitProfileMinuit(yFitIC, wsJoYAvg, wsResSum)
     fitProfileMantidFit(yFitIC, wsJoYAvg, wsResSum)
 
-    printYSpaceFitResults(wsJoYAvg.name())
+    printYSpaceFitResults()
 
     yfitResults = ResultsYFitObject(IC, yFitIC, wsTOF.name(), wsJoYAvg.name())
     yfitResults.save()
@@ -104,7 +107,10 @@ def subtract_profiles_except_lightest(ic, ws):
     if len(ic.masses) == 1:
         return
 
-    ws_name_lightest_mass = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_' + str(min(ic.masses)) + '_ncp'
+    # TODO: Make the fetching of these workspaces more robust, prone to error
+    profiles_table = mtd[ic.name + '_initial_parameters']
+    lightest_mass_str = profiles_table.column('label')[np.argmin(profiles_table.column('mass'))]
+    ws_name_lightest_mass = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_' + lightest_mass_str + '_ncp'
     ws_name_profiles = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_total_ncp'
 
     wsNcpExceptFirst = Minus(mtd[ws_name_profiles], mtd[ws_name_lightest_mass], 
@@ -127,7 +133,9 @@ def switchFirstTwoAxis(A):
 def ySpaceReduction(wsTOF, mass0, yFitIC, ic):
     """Seperate procedures depending on masking specified."""
 
-    ws_name_lightest_mass = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_' + str(min(ic.masses)) + '_ncp'
+    profiles_table = mtd[ic.name + '_initial_parameters']
+    lightest_mass_str = profiles_table.column('label')[np.argmin(profiles_table.column('mass'))]
+    ws_name_lightest_mass = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_' + lightest_mass_str + '_ncp'
     ncp = mtd[ws_name_lightest_mass].extractY()
 
     rebinPars = yFitIC.range_for_rebinning_in_y_space
@@ -199,7 +207,7 @@ def replaceZerosWithNCP(ws, ncp):
     ]  # mask of ncp adjusted for last col present or not
 
     wsMasked = CloneWorkspace(ws, OutputWorkspace=ws.name() + "_NCPMasked")
-    passDataIntoWS(dataX, dataY, dataE, wsMasked)
+    pass_data_into_ws(dataX, dataY, dataE, wsMasked)
     SumSpectra(wsMasked, OutputWorkspace=wsMasked.name() + "_Sum")
     return wsMasked
 
@@ -249,7 +257,7 @@ def dataXBining(ws, xp):
     dataE[dataY == 0] = 0
 
     wsXBins = CloneWorkspace(ws, OutputWorkspace=ws.name() + "_XBinned")
-    wsXBins = passDataIntoWS(dataX, dataY, dataE, wsXBins)
+    wsXBins = pass_data_into_ws(dataX, dataY, dataE, wsXBins)
     return wsXBins
 
 
@@ -414,15 +422,6 @@ def extractWS(ws):
     return ws.extractX(), ws.extractY(), ws.extractE()
 
 
-def passDataIntoWS(dataX, dataY, dataE, ws):
-    "Modifies ws data to input data"
-    for i in range(ws.getNumberHistograms()):
-        ws.dataX(i)[:] = dataX[i, :]
-        ws.dataY(i)[:] = dataY[i, :]
-        ws.dataE(i)[:] = dataE[i, :]
-    return ws
-
-
 def symmetrizeWs(avgYSpace):
     """
     Symmetrizes workspace after weighted average.
@@ -438,7 +437,7 @@ def symmetrizeWs(avgYSpace):
         dataYS, dataES = weightedSymArr(dataY, dataE)
 
     wsSym = CloneWorkspace(avgYSpace, OutputWorkspace=avgYSpace.name() + "_sym")
-    wsSym = passDataIntoWS(dataX, dataYS, dataES, wsSym)
+    wsSym = pass_data_into_ws(dataX, dataYS, dataES, wsSym)
     return wsSym
 
 
@@ -800,9 +799,9 @@ def model(x, A, sig_x, sig_y, sig_z):
             "Fitting Model not recognized, available options: 'gauss', 'gcc4c6', 'gcc4', 'gcc6', 'ansiogauss' gauss3d'"
         )
 
-    print("\nShared Parameters: ", [key for key in sharedPars])
-    print(
-        "\nUnshared Parameters: ", [key for key in defaultPars if key not in sharedPars]
+    logger.notice(f"\nShared Parameters: {[key for key in sharedPars]}")
+    logger.notice(
+        f"\nUnshared Parameters: {[key for key in defaultPars if key not in sharedPars]}"
     )
 
     assert all(
@@ -903,6 +902,8 @@ def createCorrelationTableWorkspace(wsYSpaceSym, parameters, corrMatrix):
         tableWS.addColumn(type="float", name=p)
     for p, arr in zip(parameters, corrMatrix):
         tableWS.addRow([p] + list(arr))
+    print_table_workspace(tableWS)
+
 
 
 def runMinos(mObj, yFitIC, constrFunc, wsName):
@@ -965,7 +966,7 @@ def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showP
     """
     # Reason for two distinct operations inside the same function is that its easier
     # to build the minos plots for each parameter as they are being calculated.
-    print("\nRunning Minos ... \n")
+    logger.notice("\nRunning Minos ... \n")
 
     # Set format of subplots
     height = 2
@@ -1224,7 +1225,7 @@ def oddPointsRes(x, res):
 
 
 def fitProfileMantidFit(yFitIC, wsYSpaceSym, wsRes):
-    print("\nFitting on the sum of spectra in the West domain ...\n")
+    logger.notice("\nFitting on the sum of spectra in the West domain ...\n")
     for minimizer in ["Levenberg-Marquardt", "Simplex"]:
         if yFitIC.fitting_model== "gauss":
             function = f"""composite=Convolution,FixResolution=true,NumDeriv=true;
@@ -1282,40 +1283,10 @@ def fitProfileMantidFit(yFitIC, wsYSpaceSym, wsRes):
     return
 
 
-def printYSpaceFitResults(wsJoYName):
-    print("\nFit in Y Space results:")
-    foundWS = []
-    try:
-        wsFitLM = mtd[wsJoYName + "_lm_Parameters"]
-        foundWS.append(wsFitLM)
-    except KeyError:
-        pass
-    try:
-        wsFitSimplex = mtd[wsJoYName + "_simplex_Parameters"]
-        foundWS.append(wsFitSimplex)
-    except KeyError:
-        pass
-    try:
-        wsFitMinuit = mtd[wsJoYName + "_minuit_Parameters"]
-        foundWS.append(wsFitMinuit)
-    except KeyError:
-        pass
-
-    for tableWS in foundWS:
-        print("\n" + " ".join(tableWS.getName().split("_")[-3:]) + ":")
-        # print("    ".join(tableWS.keys()))
-        for key in tableWS.keys():
-            if key == "Name":
-                print(
-                    f"{key:>20s}:  "
-                    + "  ".join([f"{elem:7.8s}" for elem in tableWS.column(key)])
-                )
-            else:
-                print(
-                    f"{key:>20s}: "
-                    + "  ".join([f"{elem:7.4f}" for elem in tableWS.column(key)])
-                )
-    print("\n")
+def printYSpaceFitResults():
+    for ws_name in mtd.getObjectNames():
+        if ws_name.endswith('Parameters'):
+            print_table_workspace(mtd[ws_name])
 
 
 class ResultsYFitObject:
@@ -1388,7 +1359,7 @@ def save(self):
 
 
 def runGlobalFit(wsYSpace, wsRes, IC, yFitIC):
-    print("\nRunning GLobal Fit ...\n")
+    logger.notice("\nRunning GLobal Fit ...\n")
 
     dataX, dataY, dataE, dataRes, instrPars = extractData(wsYSpace, wsRes, IC)
     dataX, dataY, dataE, dataRes, instrPars = takeOutMaskedSpectra(
@@ -1418,7 +1389,7 @@ def runGlobalFit(wsYSpace, wsRes, IC, yFitIC):
     # Minuit Fit with global cost function and local+global parameters
     initPars = minuitInitialParameters(defaultPars, sharedPars, len(dataY))
 
-    print("\nRunning Global Fit ...\n")
+    logger.notice("\nRunning Global Fit ...\n")
     m = Minuit(totCost, **initPars)
 
     for i in range(len(dataY)):  # Set limits for unshared parameters
@@ -1467,7 +1438,7 @@ def constr(*pars):
         m.scipy(constraints=optimize.NonlinearConstraint(constr, 0, np.inf))
 
     t1 = time.time()
-    print(f"\nTime of fitting: {t1-t0:.2f} seconds")
+    logger.notice(f"\nTime of fitting: {t1-t0:.2f} seconds")
 
     # Explicitly calculate errors
     m.hesse()
@@ -1528,7 +1499,7 @@ def groupDetectors(ipData, yFitIC):
 
     checkNGroupsValid(yFitIC, ipData)
 
-    print(f"\nNumber of gropus: {yFitIC.number_of_global_fit_groups}")
+    logger.notice(f"\nNumber of gropus: {yFitIC.number_of_global_fit_groups}")
 
     L1 = ipData[:, -1].copy()
     theta = ipData[:, 2].copy()
@@ -1661,11 +1632,11 @@ def formIdxList(clusters):
         idxList.append(list(idxs))
 
     # Print groupings information
-    print("\nGroups formed successfully:\n")
+    logger.notice("\nGroups formed successfully:\n")
     groupLen = np.array([len(group) for group in idxList])
     unique, counts = np.unique(groupLen, return_counts=True)
     for length, no in zip(unique, counts):
-        print(f"{no} groups with {length} detectors.")
+        logger.notice(f"{no} groups with {length} detectors.")
 
     return idxList
 
@@ -1856,19 +1827,18 @@ def create_table_for_global_fit_parameters(wsName, m, chi2):
     t.addColumn(type="float", name="Value")
     t.addColumn(type="float", name="Error")
 
-    print(f"Value of Chi2/ndof: {chi2:.2f}")
-    print(f"Migrad Minimum valid: {m.valid}")
-    print("\nResults of Global Fit:\n")
+    logger.notice(f"Value of Chi2/ndof: {chi2:.2f}")
+    logger.notice(f"Migrad Minimum valid: {m.valid}")
     for p, v, e in zip(m.parameters, m.values, m.errors):
-        print(f"{p:>7s} = {v:>8.4f} \u00B1 {e:<8.4f}")
         t.addRow([p, v, e])
 
     t.addRow(["Cost function", chi2, 0])
+    print_table_workspace(t)
 
 
 def plotGlobalFit(dataX, dataY, dataE, mObj, totCost, wsName, yFitIC):
     if len(dataY) > 10:
-        print("\nToo many axes to show in figure, skipping the plot ...\n")
+        logger.notice("\nToo many axes to show in figure, skipping the plot ...\n")
         return
 
     rows = 2
@@ -1910,5 +1880,5 @@ def plotGlobalFit(dataX, dataY, dataE, mObj, totCost, wsName, yFitIC):
 
 def save_workspaces(yFitIC):
     for ws_name in mtd.getObjectNames():
-        if ws_name.endswith('Parameters') or ws_name.endswith('Workspace'):
+        if ws_name.endswith('Parameters') or ws_name.endswith('parameters') or ws_name.endswith('Workspace'):
             SaveAscii(ws_name, str(yFitIC.figSavePath.parent / "output_files" / ws_name))

src/mvesuvio/analysis_reduction.py

@@ -13,7 +13,7 @@
                             CreateWorkspace
 
 from mvesuvio.util.analysis_helpers import numerical_third_derivative, load_resolution, load_instrument_params, \
-                                            extend_range_of_array
+                                            extend_range_of_array, print_table_workspace
 
 np.set_printoptions(suppress=True, precision=4, linewidth=200)
 
@@ -238,13 +238,13 @@ def _initialize_table_fit_parameters(self):
             OutputWorkspace=self._workspace_being_fit.name()+ "_fit_results"
         )
         table.setTitle("SciPy Fit Parameters")
-        table.addColumn(type="float", name="Spectrum")
+        table.addColumn(type="float", name="Spec")
         for label in self._profiles_table.column("label"):
-            table.addColumn(type="float", name=f"{label} intensity")
-            table.addColumn(type="float", name=f"{label} width")
-            table.addColumn(type="float", name=f"{label} center ")
-        table.addColumn(type="float", name="normalised chi2")
-        table.addColumn(type="float", name="no of iterations")
+            table.addColumn(type="float", name=f"{label} i")
+            table.addColumn(type="float", name=f"{label} w")
+            table.addColumn(type="float", name=f"{label} c")
+        table.addColumn(type="float", name="NChi2")
+        table.addColumn(type="float", name="Iter")
         return table
 
 
@@ -327,6 +327,7 @@ def _fit_neutron_compton_profiles(self):
             self._row_being_fit += 1
 
         assert np.any(self._fit_parameters), "Fitting parameters cannot be zero for all spectra!"
+        print_table_workspace(self._table_fit_results)
         return
 
 
@@ -420,8 +421,8 @@ def _set_means_and_std(self):
         intensities = np.zeros(widths.shape)
 
         for i, label in enumerate(self._profiles_table.column("label")):
-            widths[i] = self._table_fit_results.column(f"{label} width")
-            intensities[i] = self._table_fit_results.column(f"{label} intensity")
+            widths[i] = self._table_fit_results.column(f"{label} w")
+            intensities[i] = self._table_fit_results.column(f"{label} i")
             self._set_means_and_std_arrays(widths, intensities)
 
         self._create_means_table()
@@ -464,7 +465,6 @@ def _create_means_table(self):
         table.addColumn(type="float", name="mean_intensity")
         table.addColumn(type="float", name="std_intensity")
 
-        self.log().notice("\nmass    mean widths    mean intensities\n")
         for label, mass, mean_width, std_width, mean_intensity, std_intensity in zip(
             self._profiles_table.column("label"),
             self._masses,
@@ -475,17 +475,18 @@ def _create_means_table(self):
         ):
             # Explicit conversion to float required to match profiles table
             table.addRow([label, float(mass), float(mean_width), float(std_width), float(mean_intensity), float(std_intensity)])
-            self.log().notice(f"{label:6s}  {mean_width:10.5f} \u00B1 {std_width:7.5f}" + \
-                f"{mean_intensity:10.5f} \u00B1 {std_intensity:7.5f}\n")
 
         self.setPropertyValue("OutputMeansTable", table.name())
+        print_table_workspace(table, precision=5)
         return table
 
 
     def _fit_neutron_compton_profiles_to_row(self):
 
         if np.all(self._dataY[self._row_being_fit] == 0):
             self._table_fit_results.addRow(np.zeros(3*self._profiles_table.rowCount()+3))
+            spectrum_number = self._instrument_params[self._row_being_fit, 0]
+            self.log().notice(f"Skip spectrum {int(spectrum_number):3d}")
             return
 
         result = scipy.optimize.minimize(
@@ -508,7 +509,7 @@ def _fit_neutron_compton_profiles_to_row(self):
         # Store results for easier access when calculating means
         self._fit_parameters[self._row_being_fit] = tableRow 
 
-        self.log().notice(' '.join(str(tableRow).split(",")).replace('[', '').replace(']', ''))
+        self.log().notice(f"Fit spectrum {int(spectrum_number):3d}: \u2713")
 
         # Pass fit profiles into workspaces
         ncp_for_each_mass, fse_for_each_mass = self._neutron_compton_profiles(fitPars)

src/mvesuvio/main/__init__.py

@@ -44,6 +44,9 @@ def __set_up_parser():
 
 
 def __setup_config(args):
+
+    __set_logging_properties()
+
     config_dir = handle_config.VESUVIO_CONFIG_PATH
     handle_config.setup_config_dir(config_dir)
     ipfolder_dir = handle_config.VESUVIO_IPFOLDER_PATH
@@ -80,10 +83,24 @@ def __setup_config(args):
     handle_config.check_dir_exists("IP folder", ipfolder_dir)
 
 
+def __set_logging_properties():
+    from mantid.kernel import ConfigService
+    ConfigService.setString("logging.loggers.root.channel.class", "SplitterChannel")
+    ConfigService.setString("logging.loggers.root.channel.channel1", "consoleChannel")
+    ConfigService.setString("logging.loggers.root.channel.channel2", "fileChannel")
+    ConfigService.setString("logging.channels.consoleChannel.class", "ConsoleChannel")
+    ConfigService.setString("logging.channels.fileChannel.class", "FileChannel")
+    ConfigService.setString("logging.channels.fileChannel.path", "mantid.log")
+    ConfigService.setString("logging.channels.fileChannel.formatter.class", "PatternFormatter")
+    ConfigService.setString("logging.channels.fileChannel.formatter.pattern", "%Y-%m-%d %H:%M:%S,%i [%I] %p %s - %t")
+    # Set properties on Mantid.user.properties not working due to Mantid bug
+    # Need to set properties on file in Mantid installation
+    mantid_properties_file = path.join(ConfigService.getPropertiesDir(), "Mantid.properties")
+    ConfigService.saveConfig(mantid_properties_file)
+    return
+
+
 def __run_analysis():
-    environ["MANTIDPROPERTIES"] = path.join(
-        handle_config.VESUVIO_CONFIG_PATH, "Mantid.user.properties"
-    )
     from mvesuvio.main.run_routine import Runner 
     Runner().run()