From 87e6f8d4d46e5d0be2125a343120300d8634c67e Mon Sep 17 00:00:00 2001
From: Ian Jackson <ijackson@chiark.greenend.org.uk>
Date: Sun, 26 Nov 2017 14:15:53 +0000
Subject: [PATCH] symbolic.py: get_python produces answer of right shape

Signed-off-by: Ian Jackson <ijackson@chiark.greenend.org.uk>
---
 helixish.py                    |  6 +++++-
 symbolic-python-lambda-test.py | 11 +++++++++++
 symbolic.py                    | 19 +++++++++++++++++++
 3 files changed, 35 insertions(+), 1 deletion(-)
 create mode 100644 symbolic-python-lambda-test.py

diff --git a/helixish.py b/helixish.py
index 3479910..4543898 100644
--- a/helixish.py
+++ b/helixish.py
@@ -7,6 +7,8 @@ from numpy import cos, sin
 import sys
 from moebdebug import dbg
 
+import symbolic
+
 def augment(v): return np.append(v, 1)
 def augment0(v): return np.append(v, 0)
 def unaugment(v): return v[0:3]
@@ -15,6 +17,8 @@ findcurve_subproc = None
 
 class HelixishCurve():
   def __init__(hc, cp):
+    symbolic.calculate()
+
     p = cp[0]
     q = cp[3]
     dp = unit_v(cp[1]-cp[0])
@@ -148,4 +152,4 @@ class HelixishCurve():
       if l is None: break
       findcurve_result = l[0:5]
 
-    
+    symbolic.get_python(something)
diff --git a/symbolic-python-lambda-test.py b/symbolic-python-lambda-test.py
new file mode 100644
index 0000000..500c4f9
--- /dev/null
+++ b/symbolic-python-lambda-test.py
@@ -0,0 +1,11 @@
+import symbolic
+
+symbolic.calculate()
+l = symbolic.get_python()
+print(repr(l))
+r = l(1,2,3,4,5,6,7)
+print(r)
+r = r[:,0]
+print(r)
+x,y,z = r
+print(x,y,z)
diff --git a/symbolic.py b/symbolic.py
index 8607034..f787079 100644
--- a/symbolic.py
+++ b/symbolic.py
@@ -4,6 +4,8 @@ import itertools
 
 from moedebug import dbg_enable
 
+from sympy.utilities.lambdify import lambdify, implemented_function
+
 r, theta, s, la, mu, kappa = symbols('r theta s lambda mu kappa')
 
 # start      original formulation
@@ -21,7 +23,12 @@ def dbg(*args):
     print('\n          =\n')
     pprint(cse(eval(vn)))
 
+calculated = False
+
 def calculate():
+  global calculated
+  if calculated: return
+
   p_start = Matrix([
     r * (1 - cos(theta)),
     r * sin(theta),
@@ -133,6 +140,7 @@ def calculate():
   global sh, th
   sh, th = symbols('alpha beta')
 
+  global q_sqparm
   q_dirn_sqparm = q_dirn_orgcoords.replace(s, sh**2).replace(t, th**2)
   q_sqparm      = q_orgcoords     .replace(s, sh**2).replace(t, th**2)
 
@@ -152,6 +160,8 @@ def calculate():
   result_dirnscaled = q_sqparm.col_join(q_dirn_dirnscaled)
   dbg('result_dirnscaled')
 
+  calculated = True
+
 params = ('sh','th','la','mu','gamma','kappa')
 
 def ourccode(*a, **kw):
@@ -231,3 +241,12 @@ def gen_C():
   gen_x_extract()
   gen_f_populate()
   gen_j_populate()
+
+def get_python():
+  # https://github.com/sympy/sympy/issues/13642
+  # "lambdify sinc gives wrong answer!"
+  sinc_fixed = Function('sinc_fixed')
+  implemented_function(sinc_fixed, lambda x: np.sinc(x/np.pi))
+  p = list(map(eval,params))
+  p.append(s)
+  return lambdify(p, q_sqparm)
-- 
2.30.2