Compare with legacy McaTheory (WIP)

PyMca5/tests/XrfTest.py

@@ -267,10 +267,8 @@ def testTrainingDataFilePresence(self):
         self.assertTrue(os.path.isfile(trainingDataFile),
                         "File %s is not an actual file" % trainingDataFile)
 
-    def testTrainingDataFit(self):
+    def _readTrainingData(self):
         from PyMca5.PyMcaIO import specfilewrapper as specfile
-        from PyMca5.PyMcaPhysics.xrf import LegacyMcaTheory
-        from PyMca5.PyMcaPhysics.xrf import ConcentrationsTool
         from PyMca5.PyMcaIO import ConfigDict
         trainingDataFile = os.path.join(self.dataDir, "XRFSpectrum.mca")
         self.assertTrue(os.path.isfile(trainingDataFile),
@@ -285,18 +283,52 @@ def testTrainingDataFit(self):
                         "Training data 1st scan should contain no MCAs")
         y = mcaData = sf[1].mca(1)
         sf = None
+        x = numpy.arange(y.size).astype(numpy.float64)
 
         # perform the actual XRF analysis
         configuration = ConfigDict.ConfigDict()
         configuration.readfp(StringIO(cfg))
+
+        return x, y, configuration
+
+    def _readStainlessSteelData(self):
+        from PyMca5.PyMcaIO import specfilewrapper as specfile
+        from PyMca5.PyMcaIO import ConfigDict
+
+        # read the data
+        dataFile = os.path.join(self.dataDir, "Steel.spe")
+        self.assertTrue(os.path.isfile(dataFile),
+                        "File %s is not an actual file" % dataFile)
+        sf = specfile.Specfile(dataFile)
+        self.assertTrue(len(sf) == 1, "File %s cannot be read" % dataFile)
+        self.assertTrue(sf[0].nbmca() == 1,
+                        "Spe file should contain MCA data")
+        y = counts = sf[0].mca(1)
+        x = channels = numpy.arange(y.size).astype(numpy.float64)
+        sf = None
+
+        # read the fit configuration
+        configFile = os.path.join(self.dataDir, "Steel.cfg")
+        self.assertTrue(os.path.isfile(configFile),
+                        "File %s is not an actual file" % configFile)
+        configuration = ConfigDict.ConfigDict()
+        configuration.read(configFile)
+        # configure the fit
+        # make sure no secondary excitations are used
+        configuration["concentrations"]["usemultilayersecondary"] = 0
+
+        return x, y, configuration
+
+    def testTrainingDataFit(self):
+        from PyMca5.PyMcaIO import specfilewrapper as specfile
+        from PyMca5.PyMcaPhysics.xrf import LegacyMcaTheory
+        from PyMca5.PyMcaPhysics.xrf import ConcentrationsTool
+
+        x, y, configuration = self._readTrainingData()
+
+        # perform the actual XRF analysis
         mcaFit = LegacyMcaTheory.LegacyMcaTheory()
-        configuration=mcaFit.configure(configuration)
-        x = numpy.arange(y.size).astype(numpy.float64)
-        mcaFit.setData(x, y,
-                       xmin=configuration["fit"]["xmin"],
-                       xmax=configuration["fit"]["xmax"])
-        mcaFit.estimate()
-        fitResult, result = mcaFit.startFit(digest=1)
+        fitResult, result = self._configAndFit(x, y, configuration, mcaFit)
 
         # fit is already done, calculate the concentrations
         concentrationsConfiguration = configuration["concentrations"]
@@ -350,39 +382,14 @@ def testTrainingDataFit(self):
                 "Error for <%s> concentration %g != %g" % (key, internal, fp))
 
     def testStainlessSteelDataFit(self):
-        from PyMca5.PyMcaIO import specfilewrapper as specfile
         from PyMca5.PyMcaPhysics.xrf import LegacyMcaTheory
         from PyMca5.PyMcaPhysics.xrf import ConcentrationsTool
-        from PyMca5.PyMcaIO import ConfigDict
 
-        # read the data
-        dataFile = os.path.join(self.dataDir, "Steel.spe")
-        self.assertTrue(os.path.isfile(dataFile),
-                        "File %s is not an actual file" % dataFile)
-        sf = specfile.Specfile(dataFile)
-        self.assertTrue(len(sf) == 1, "File %s cannot be read" % dataFile)
-        self.assertTrue(sf[0].nbmca() == 1,
-                        "Spe file should contain MCA data")
-        y = counts = sf[0].mca(1)
-        x = channels = numpy.arange(y.size).astype(numpy.float64)
-        sf = None
+        x, y, configuration = self._readStainlessSteelData()
 
-        # read the fit configuration
-        configFile = os.path.join(self.dataDir, "Steel.cfg")
-        self.assertTrue(os.path.isfile(configFile),
-                        "File %s is not an actual file" % configFile)
-        configuration = ConfigDict.ConfigDict()
-        configuration.read(configFile)
         # configure the fit
-        # make sure no secondary excitations are used
-        configuration["concentrations"]["usemultilayersecondary"] = 0
         mcaFit = LegacyMcaTheory.LegacyMcaTheory()
-        configuration=mcaFit.configure(configuration)
-        mcaFit.setData(x, y,
-                       xmin=configuration["fit"]["xmin"],
-                       xmax=configuration["fit"]["xmax"])
-        mcaFit.estimate()
-        fitResult, result = mcaFit.startFit(digest=1)
+        fitResult, result = self._configAndFit(x, y, configuration, mcaFit)
 
         # concentrations
         # fit is already done, calculate the concentrations
@@ -456,12 +463,7 @@ def testStainlessSteelDataFit(self):
         # in order to get the good fundamental parameters
         configuration["concentrations"]['usematrix'] = 1
         configuration["concentrations"]["usemultilayersecondary"] = 2
-        mcaFit.setConfiguration(configuration)
-        mcaFit.setData(x, y,
-                       xmin=configuration["fit"]["xmin"],
-                       xmax=configuration["fit"]["xmax"])
-        mcaFit.estimate()
-        fitResult, result = mcaFit.startFit(digest=1)
+        fitResult, result = self._configAndFit(x, y, configuration, mcaFit)
 
         # concentrations
         # fit is already done, calculate the concentrations
@@ -511,12 +513,7 @@ def testStainlessSteelDataFit(self):
                                                          "Ni", "-", "-",
                                                          "-","-","-"]
         mcaFit = LegacyMcaTheory.LegacyMcaTheory()
-        configuration=mcaFit.configure(configuration)
-        mcaFit.setData(x, y,
-                       xmin=configuration["fit"]["xmin"],
-                       xmax=configuration["fit"]["xmax"])
-        mcaFit.estimate()
-        fitResult, result = mcaFit.startFit(digest=1)
+        fitResult, result = self._configAndFit(x, y, configuration, mcaFit)
 
         # concentrations
         # fit is already done, calculate the concentrations
@@ -546,6 +543,79 @@ def testStainlessSteelDataFit(self):
                 "Strategy: Element %s discrepancy too large %.1f %%" % \
                   (element.split()[0], delta))
 
+    def testCompareLegacyMcaTheory(self):
+        x, y, configuration = self._readTrainingData()
+        self._testCompareLegacyMcaTheory(x, y, configuration)
+        return
+        x, y, configuration = self._readStainlessSteelData()
+        self._testCompareLegacyMcaTheory(x, y, configuration)
+
+    def _testCompareLegacyMcaTheory(self, x, y, configuration):
+        from PyMca5.PyMcaPhysics.xrf import LegacyMcaTheory
+        from PyMca5.PyMcaPhysics.xrf import NewClassMcaTheory
+
+        mcaFitLegacy = LegacyMcaTheory.LegacyMcaTheory()
+        mcaFit = NewClassMcaTheory.McaTheory()
+        fitResult1, result1 = self._configAndFit(
+            x, y, configuration, mcaFitLegacy, tmpflag=True)
+        fitResult2, result2 = self._configAndFit(
+            x, y, configuration, mcaFit, tmpflag=True)
+
+        # Compare internals
+        self.assertEqual(mcaFitLegacy.config, mcaFit.config)
+
+        self.assertEqual(mcaFitLegacy._fluoRates, mcaFit._fluoRates)
+
+        # Compare line groups
+        linegroups1 = mcaFitLegacy.PEAKS0
+        linegroups2 = mcaFit._lineGroups
+        linegroups1 = [[[line[1], line[0], name]
+                        for name, line in zip(names, lines)] 
+                        for names, lines in zip(mcaFitLegacy.PEAKS0NAMES, linegroups1)]
+        nsource = mcaFit._nSourceLines
+        if nsource:
+            linegroups2 = linegroups2[:-nsource]
+        # self.assertEqual(linegroups1, linegroups2)
+        self.assertEqual(len(linegroups1), len(linegroups2))
+        for lg1, lg2 in zip(linegroups1, linegroups2):
+            self.assertEqual(len(lg1), len(lg2))
+            for line1, line2 in zip(lg1, lg2):
+                self.assertEqual(line1[0], line2[0])
+                numpy.testing.assert_allclose(line1[1], line2[1], rtol=1e-9)
+                self.assertEqual(line1[2], line2[2])
+
+        # Compare escape line groups
+        linegroups1 = mcaFitLegacy.PEAKS0ESCAPE
+        linegroups2 = mcaFit._escapeLineGroups
+        nsource = mcaFit._nSourceLines
+        if nsource:
+            linegroups2 = linegroups2[:-nsource]
+        #self.assertEqual(linegroups1, linegroups2)
+        self.assertEqual(len(linegroups1), len(linegroups2))
+        for lg1, lg2 in zip(linegroups1, linegroups2):
+            self.assertEqual(len(lg1), len(lg2))
+            for elines1, elines2 in zip(lg1, lg2):
+                self.assertEqual(len(elines1), len(elines2))
+                for line1, line2 in zip(elines1, elines2):
+                    self.assertEqual(line1[0], line2[0])
+                    numpy.testing.assert_allclose(line1[1], line2[1], rtol=1e-9)
+                    self.assertEqual(line1[2], line2[2])
+
+        # Compare fit results
+        self.assertEqual(fitResult1, fitResult2)
+        self.assertEqual(result1, result2)
+
+    def _configAndFit(self, x, y, configuration, mcaFit, tmpflag=False):
+        configuration = mcaFit.configure(configuration)
+        if tmpflag:
+            return None, None
+        mcaFit.setData(x, y,
+                       xmin=configuration["fit"]["xmin"],
+                       xmax=configuration["fit"]["xmax"])
+        mcaFit.estimate()
+        return mcaFit.startFit(digest=1)
+
+
 def getSuite(auto=True):
     testSuite = unittest.TestSuite()
     if auto: