From df35cb174bbce3bda4203365909d5c6829649e59 Mon Sep 17 00:00:00 2001
From: Stephen Mardson McQuay <stephen@mcquay.me>
Date: Wed, 5 May 2010 23:03:12 -0600
Subject: [PATCH] made a get_error function that is generic and functional. I'm
 still getting unsavory numbers for the cubic, 2-D case on my single test
 case. I will have to try out more refined grids, albeit the quadratic is
 consistently improving results for my case.

---
 bin/driver.py              |  2 +-
 bin/test2d-connectivity.py |  3 +-
 bin/test_baker.py          |  7 ++---
 lib/baker/__init__.py      | 59 ++++++++++++++++++++++++++++++++++----
 4 files changed, 59 insertions(+), 12 deletions(-)

diff --git a/bin/driver.py b/bin/driver.py
index 6d2a35a..db5548a 100755
--- a/bin/driver.py
+++ b/bin/driver.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python
 
 import sys
 from optparse import OptionParser
diff --git a/bin/test2d-connectivity.py b/bin/test2d-connectivity.py
index 524376d..7e2b818 100755
--- a/bin/test2d-connectivity.py
+++ b/bin/test2d-connectivity.py
@@ -37,7 +37,7 @@ if __name__ == '__main__':
     sys.exit(1)
 
 
-  print "total points", total_points
+  smblog.debug("total points: %d" % total_points)
 
 
 
@@ -67,3 +67,4 @@ if __name__ == '__main__':
     good.append(d[True])
 
   draw_gb(bad = bad, good = good)
+
diff --git a/bin/test_baker.py b/bin/test_baker.py
index a6d1eb4..0bff1f2 100755
--- a/bin/test_baker.py
+++ b/bin/test_baker.py
@@ -1,6 +1,5 @@
 #!/usr/bin/env python
 
-import math
 from baker        import run_baker
 from baker.tools import improved_answer
 
@@ -9,10 +8,12 @@ from baker.tools import smblog
 from grid.DD      import grid
 from grid.simplex import contains
 
+import numpy as np
+
 def exact_func(X):
   x = X[0]
   y = X[0]
-  return 1 - math.sin((x-0.5)**2 + (y-0.5)**2)
+  return 1 - np.sin(np.power((x-0.5), 2) + np.power((y-0.5), 2))
 
 if __name__ == '__main__':
   points = [
@@ -33,8 +34,6 @@ if __name__ == '__main__':
   g.construct_connectivity()
   R = g.create_mesh(range(3))
 
-  print contains(X, R.points)
-
   try:
     extra = int(sys.argv[1])
   except:
diff --git a/lib/baker/__init__.py b/lib/baker/__init__.py
index c16522d..323cf68 100644
--- a/lib/baker/__init__.py
+++ b/lib/baker/__init__.py
@@ -88,7 +88,7 @@ def qlinear(X, R):
   """
 
   phis = get_phis(X, R.points)
-  qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)])
+  qlin = np.sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)])
   return phis, qlin
 
 def qlinear_3D(X, R):
@@ -185,6 +185,54 @@ def get_error_cubic(phi, R, S):
 
   return error_term
 
+def get_error_sauron(phi, R, S, order = 2):
+  smblog.debug("len(phi): %d"% len(phi))
+  B = [] # baker eq 9
+  w = [] # baker eq 11
+
+  p  = pattern(order, len(phi), offset = -1)
+  smblog.debug("pattern: %s" % p)
+
+  for (s,q) in zip(S.points, S.q):
+    cur_phi, cur_qlin = qlinear(s, R)
+    l = []
+    for i in p:
+      cur_sum = cur_phi[i[0]]
+      for j in i[1:]:
+        cur_sum *= cur_phi[j]
+      l.append(cur_sum)
+
+    B.append(l)
+    w.append(q - cur_qlin)
+
+  smblog.debug("B: %s" % B)
+  smblog.debug("w: %s" % w)
+    
+
+  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:
+    abc = np.linalg.solve(A,b)
+  except np.linalg.LinAlgError as e:
+    smblog.error("linear calculation went bad, resorting to np.linalg.pinv: %s" % e)
+    abc = np.dot(np.linalg.pinv(A), b)
+
+  smblog.debug(len(abc) == len(p))
+
+  error_term = 0.0
+  for (a, i) in zip(abc, p):
+    cur_sum = a
+    for j in i:
+      cur_sum *= phi[j]
+    error_term += cur_sum
+    
+  smblog.debug("error_term smb: %s" % error_term)
+  return error_term, abc
 
 def run_baker(X, R, S, order=2):
   """
@@ -196,6 +244,7 @@ def run_baker(X, R, S, order=2):
     R = Simplex
     S = extra points
   """
+  smblog.debug("order = %d" % order)
 
   answer = {
              'qlin': None,
@@ -208,12 +257,10 @@ def run_baker(X, R, S, order=2):
   if order == 1:
     answer['qlin'] = qlin
     return answer
-  elif order == 2:
-    error_term = get_error_quadratic(phi, R, S)
-  elif order == 3:
-    error_term = get_error_cubic(phi, R, S)
+  elif order in (2,3):
+    error_term, abc = get_error_sauron(phi, R, S, order)
   else:
-    raise smberror('unacceptable order for baker method')
+    raise smberror('unsupported order for baker method')
 
   q_final = qlin + error_term