moved the quadratic error calculation out of the baker method, and added a cubic version too

lib/baker/__init__.py

@@ -102,6 +102,88 @@ def qlinear_3D(X, R):
   qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)])
   return phis, qlin
 
+def get_error_quadratic(phi, R, S):
+  B = [] # baker eq 9
+  w = [] # baker eq 11
+
+  for (s, q) in zip(S.points, S.q):
+    cur_phi, cur_qlin = qlinear(s, R)
+    (phi1, phi2, phi3) = cur_phi
+
+    B.append(
+              [
+                phi1 * phi2,
+                phi2 * phi3,
+                phi3 * phi1,
+              ]
+            )
+    w.append(q - cur_qlin)
+
+  B = np.array(B)
+  w = np.array(w)
+
+  A = np.dot(B.T, B)
+  b = np.dot(B.T, w)
+
+  # baker solve eq 10
+  try:
+    (a, b, c) = np.linalg.solve(A,b)
+  except np.linalg.LinAlgError as e:
+    print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv", e
+    (a, b, c) = np.dot(np.linalg.pinv(A), b)
+
+  error_term =  a * phi[0] * phi[1]\
+              + b * phi[1] * phi[2]\
+              + c * phi[2] * phi[0]
+
+  return error_term, a, b, c
+
+def get_error_cubic(phi, R, S):
+  B = [] # baker eq 9
+  w = [] # baker eq 11
+
+  for (s, q) in zip(S.points, S.q):
+    cur_phi, cur_qlin = qlinear(s, R)
+    (phi1, phi2, phi3) = cur_phi
+
+    # basing this on eq 17
+    B.append(
+              [
+                phi1 * phi2,        # a
+                phi1 * phi3,        # b
+                phi2 * phi1,        # c
+                phi2 * phi3,        # d
+                phi3 * phi1,        # e
+                phi3 * phi2,        # f
+                phi1 * phi2 * phi3, # g
+              ]
+            )
+    w.append(q - cur_qlin)
+
+  B = np.array(B)
+  w = np.array(w)
+
+  A = np.dot(B.T, B)
+  b = np.dot(B.T, w)
+
+  # baker solve eq 10
+  try:
+    (a, b, c, d, e, f, g) = np.linalg.solve(A,b)
+  except np.linalg.LinAlgError as e:
+    print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv", e
+    (a, b, c, d, e, f, g) = np.dot(np.linalg.pinv(A), b)
+
+  error_term =  a * phi[0] * phi[1]\
+              + b * phi[0] * phi[2]\
+              + c * phi[1] * phi[0]\
+              + d * phi[1] * phi[2]\
+              + e * phi[2] * phi[0]\
+              + f * phi[2] * phi[1]\
+              + g * phi[0] * phi[1] * phi[2]\
+
+  return error_term, a, b, c
+
+
 def run_baker(X, R, S, order=2):
   """
     This is the main function to call to get an interpolation to X from the input meshes
@@ -135,38 +217,12 @@ def run_baker(X, R, S, order=2):
               }
     return answer
 
-  B = [] # baker eq 9
-  w = [] # baker eq 11
-
-  for (s, q) in zip(S.points, S.q):
-    cur_phi, cur_qlin = qlinear(s, R)
-    (phi1, phi2, phi3) = cur_phi
-
-    B.append(
-              [
-                phi1 * phi2,
-                phi2 * phi3,
-                phi3 * phi1,
-              ]
-            )
-    w.append(q - cur_qlin)
-
-  B = np.array(B)
-  w = np.array(w)
-
-  A = np.dot(B.T, B)
-  b = np.dot(B.T, w)
-
-  # baker solve eq 10
-  try:
-    (a, b, c) = np.linalg.solve(A,b)
-  except np.linalg.LinAlgError as e:
-    print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv", e
-    (a, b, c) = np.dot(np.linalg.pinv(A), b)
-
-  error_term =  a * phi[0] * phi[1]\
-              + b * phi[1] * phi[2]\
-              + c * phi[2] * phi[0]
+  if order == 2:
+    error_term, a, b, c = get_error_quadratic(phi, R, S)
+  elif order == 3:
+    error_term, a, b, c = get_error_cubic(phi, R, S)
+  else:
+    raise smberror('unacceptable order for baker method')
 
   q_final = qlin + error_term