Merge pull request #24 from rchristie/calc_error

Calculate RMS and max projection errors

src/scaffoldfitter/fitter.py

@@ -682,6 +682,27 @@ def getDataWeightField(self):
     def getActiveDataNodesetGroup(self):
         return self._activeDataNodesetGroup
 
+    def getDataRMSAndMaximumProjectionError(self):
+        """
+        Get RMS and maximum error for all active data and marker point projections.
+        No group weights are applied.
+        :return: RMS error value, maximum error value. Values are None if there is no data or bad fields.
+        """
+        with ChangeManager(self._fieldmodule):
+            error = self._fieldmodule.createFieldMagnitude(self._dataDeltaField)
+            msError = self._fieldmodule.createFieldNodesetMeanSquares(error, self._activeDataNodesetGroup)
+            rmsError = self._fieldmodule.createFieldSqrt(msError)
+            maxError = self._fieldmodule.createFieldNodesetMaximum(error, self._activeDataNodesetGroup)
+            fieldcache = self._fieldmodule.createFieldcache()
+            rmsResult, rmsErrorValue = rmsError.evaluateReal(fieldcache, 1)
+            maxResult, maxErrorValue = maxError.evaluateReal(fieldcache, 1)
+            del fieldcache
+            del maxError
+            del rmsError
+            del msError
+            del error
+        return rmsErrorValue if (rmsResult == RESULT_OK) else None, maxErrorValue if (maxResult == RESULT_OK) else None
+
     def getMarkerDataFields(self):
         """
         Only call if markerGroup exists.
@@ -1077,6 +1098,10 @@ def calculateDataProjections(self, fitterStep: FitterStep):
                           " data points with data coordinates have not been projected")
                 del unprojectedDatapoints
 
+                # Write projection error measures
+                rmsErrorValue, maxErrorValue = self.getDataRMSAndMaximumProjectionError()
+                print("Data projection RMS error", rmsErrorValue, "Max error", maxErrorValue)
+
             # remove temporary objects before ChangeManager exits
             del fieldcache
 

tests/resources/square.exf

@@ -0,0 +1,109 @@
+EX Version: 3
+Region: /
+!#nodeset nodes
+Define node template: node1
+Shape. Dimension=0
+#Fields=1
+1) coordinates, coordinate, rectangular cartesian, real, #Components=3
+ x. #Values=1 (value)
+ y. #Values=1 (value)
+ z. #Values=1 (value)
+Node template: node1
+Node: 1
+ 0.000000000000000e+00
+ 0.000000000000000e+00
+ 0.000000000000000e+00
+Node: 2
+ 1.000000000000000e+00
+ 0.000000000000000e+00
+ 0.000000000000000e+00
+Node: 3
+ 0.000000000000000e+00
+ 1.000000000000000e+00
+ 0.000000000000000e+00
+Node: 4
+ 1.000000000000000e+00
+ 1.000000000000000e+00
+ 0.000000000000000e+00
+Define node template: node2
+Shape. Dimension=0
+#Fields=2
+1) marker_location, field, rectangular cartesian, element_xi, #Components=1, host mesh=mesh2d, host mesh dimension=2
+ 1. #Values=1 (value)
+2) marker_name, field, rectangular cartesian, string, #Components=1
+ 1. #Values=1 (value)
+Node template: node2
+Node: 5
+ 1  5.000000000000000e-01  5.000000000000000e-01
+ middle
+!#mesh mesh1d, dimension=1, nodeset=nodes
+Define element template: element1
+Shape. Dimension=1, line
+#Scale factor sets=0
+#Nodes=0
+#Fields=0
+Element template: element1
+Element: 1
+Element: 2
+Element: 3
+Element: 4
+Element: 5
+Element: 6
+Element: 7
+!#mesh mesh2d, dimension=2, face mesh=mesh1d, nodeset=nodes
+Define element template: element2
+Shape. Dimension=2, line*line
+#Scale factor sets=0
+#Nodes=4
+#Fields=1
+1) coordinates, coordinate, rectangular cartesian, real, #Components=3
+ x. l.Lagrange*l.Lagrange, no modify, standard node based.
+  #Nodes=4
+  1. #Values=1
+   Value labels: value
+  2. #Values=1
+   Value labels: value
+  3. #Values=1
+   Value labels: value
+  4. #Values=1
+   Value labels: value
+ y. l.Lagrange*l.Lagrange, no modify, standard node based.
+  #Nodes=4
+  1. #Values=1
+   Value labels: value
+  2. #Values=1
+   Value labels: value
+  3. #Values=1
+   Value labels: value
+  4. #Values=1
+   Value labels: value
+ z. l.Lagrange*l.Lagrange, no modify, standard node based.
+  #Nodes=4
+  1. #Values=1
+   Value labels: value
+  2. #Values=1
+   Value labels: value
+  3. #Values=1
+   Value labels: value
+  4. #Values=1
+   Value labels: value
+Element template: element2
+Element: 1
+ Faces:
+ 1 2 3 4
+ Nodes:
+ 1 2 3 4
+Group name: square
+!#nodeset nodes
+Node group:
+1..4
+!#mesh mesh1d, dimension=1, nodeset=nodes
+Element group:
+1..4
+!#mesh mesh2d, dimension=2, face mesh=mesh1d, nodeset=nodes
+Element group:
+1
+Group name: marker
+!#nodeset nodes
+Node group:
+5

tests/resources/square_error_data.exf

@@ -0,0 +1,46 @@
+EX Version: 3
+Region: /
+!#nodeset datapoints
+Define node template: node1
+Shape. Dimension=0
+#Fields=1
+1) data_coordinates, coordinate, rectangular cartesian, real, #Components=3
+ x. #Values=1 (value)
+ y. #Values=1 (value)
+ z. #Values=1 (value)
+Node template: node1
+Node: 1
+  2.000000000000000e-01
+  2.000000000000000e-01
+  5.000000000000000e-01
+Node: 2
+  8.000000000000000e-01
+  7.000000000000000e-01
+  3.000000000000000e-01
+Node: 3
+  1.100000000000000e+00
+  9.000000000000000e-01
+  1.000000000000000e-01
+Define node template: node2
+Shape. Dimension=0
+#Fields=2
+1) marker_data_coordinates, coordinate, rectangular cartesian, real, #Components=3
+ x. #Values=1 (value)
+ y. #Values=1 (value)
+ z. #Values=1 (value)
+2) marker_data_name, field, rectangular cartesian, string, #Components=1
+ 1. #Values=1 (value)
+Node template: node2
+Node: 4
+  7.000000000000000e-01
+  3.000000000000000e-01
+ -2.000000000000000e-01
+ middle
+Group name: marker
+!#nodeset datapoints
+Node group:
+4
+Group name: square
+!#nodeset datapoints
+Node group:
+1..3

tests/test_fit2d.py

@@ -57,6 +57,20 @@ def test_fit_breast2d(self):
         self.assertEqual(result, RESULT_OK)
         self.assertAlmostEqual(surfaceArea, 104501.36293993103, delta=1.0E-1)
 
+    def test_projection_error(self):
+        """
+        Test data projection RMS and maximum error calculations.
+        """
+        zinc_model_file = os.path.join(here, "resources", "square.exf")
+        zinc_data_file = os.path.join(here, "resources", "square_error_data.exf")
+        fitter = Fitter(zinc_model_file, zinc_data_file)
+        fitter.setDiagnosticLevel(1)
+        fitter.load()
+        rmsErrorValue, maxErrorValue = fitter.getDataRMSAndMaximumProjectionError()
+        TOL = 1.0E-10
+        self.assertAlmostEqual(rmsErrorValue, 0.34641016151377546, delta=TOL)  # sqrt(0.12)
+        self.assertAlmostEqual(maxErrorValue, 0.5, delta=TOL)
+
 
 if __name__ == "__main__":
     unittest.main()