X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?a=blobdiff_plain;f=symbolic.py;h=8680378ea06d906cbf58b4d230f09451e5ec281d;hb=c1001cc933c7a48935777908ba8338673cd091dc;hp=06ef2c6937f485603183389036b2ae73a0450301;hpb=66615c05daa45c98528cc6072dddc8e1cd8c1754;p=moebius3.git

diff --git a/symbolic.py b/symbolic.py
index 06ef2c6..8680378 100644
--- a/symbolic.py
+++ b/symbolic.py
@@ -9,6 +9,10 @@ import numpy as np
 
 from moedebug import *
 
+# When cse() is called for something containing BaseVector, it
+# produces infinite recursion.
+#def cse(x, *a, **kw): return ((), (x,))
+
 def dprint(*args):
   if not dbg_enabled(): return
   print(*args)
@@ -21,6 +25,8 @@ def dbg(*args):
     print('\n          =\n')
     pprint(cse(eval(vn)))
 
+def sqnorm(v): return v & v
+
 N = CoordSysCartesian('N')
 
 calculated = False
@@ -38,42 +44,105 @@ A, B, C, D = vector_symbols('A B C D')
 p = vector_symbols('p')
 
 E, H = vector_symbols('E H')
+F0, G0 = vector_symbols('F0 G0')
 En, Hn = vector_symbols('En Hn')
 
-EFl, GHl = symbols('EFl GHl')
+EFlq, HGlq = symbols('EFlq HGlq')
 
 def vector_component(v, ix):
   return v.components[N.base_vectors()[ix]]
 
 # x array in numerical algorithm has:
 #    N x 3    coordinates of points 0..N-3
-#    1        EFl = length parameter |EF| for point 1
-#    1        GHl = length parameter |GH| for point N-2
+#    1        EFlq = sqrt of length parameter |EF| for point 1
+#    1        HGlq = sqrt of length parameter |HG| for point N-2
 
 # fixed array in numerical algorithm has:
 #    4 x 3    E, En, H, Hn
 
 #def subst_vect():
 
+iterations = []
+
+class SomeIteration():
+  def __init__(ar, names, size, expr):
+    ar.names_string = names
+    ar.names = names.split(' ')
+    ar.name = ar.names[0]
+    ar.size = size
+    ar.expr = expr
+    if dbg_enabled():
+      print('\n   ' + ar.name + '\n')
+      print(expr)
+    iterations.append(ar)
+
+  def gen_calculate_cost(ar):
+    ar._gen_array()
+    cprint('for (P=0; P<(%s); P++) {' % ar.size)
+    ar._cassign()
+    cprint('}')
+
+class ScalarArray(SomeIteration):
+  def _gen_array(ar):
+    cprint('double A_%s[%s];' % (ar.name, ar.size))
+  def gen_references(ar):
+    for ai in range(0, len(ar.names)):
+      ar._gen_reference(ai, ar.names[ai])
+  def _gen_reference(ar, ai, an):
+    cprintraw('#define %s A_%s[P%+d]' % (an, ar.name, ai))
+  def _cassign(ar):
+    cassign(ar.expr, ar.name, 'tmp_'+ar.name)
+  def s(ar):
+    return symbols(ar.names_string)
+
+class CoordArray(ScalarArray):
+  def _gen_array(ar):
+    cprint('double A_%s[%s][3];' % (ar.name, ar.size))
+  def _gen_reference(ar, ai, an):
+    ScalarArray._gen_reference(ar, ai, an)
+    gen_point_coords_macro(an)
+  def _cassign(ar):
+    cassign_vector(ar.expr, ar.name, 'tmp_'+ar.name)
+  def s(ar):
+    return vector_symbols(ar.names_string)
+
+class CostComponent(SomeIteration):
+  def __init__(cc, size, expr):
+    cc.size = size
+    cc.expr = expr
+    iterations.append(cc)
+  def gen_references(cc): pass
+  def _gen_array(cc): pass
+  def _cassign(cc):
+    cassign(cc.expr, 'P_cost', 'tmp_cost')
+    cprint('cost += P_cost;')
+
 def calculate():
   global calculated
   if calculated: return
 
-  Q = B + 0.5 * (B - A)
-  R = C + 0.5 * (C - D)
-  QR = R - Q
-  BC = C - B
-  cost_ABCD = (QR & QR) / (BC & BC)
-  global cost_ABCD
-  dbg('cost_ABCD')
-  dprint(A)
+  # ---------- actual cost computation formulae ----------
+
+  global F, G
+  F = E + En * pow(EFlq, 2)
+  G = H + Hn * pow(HGlq, 2)
+
+  global a,b, al,bl, au,bu
+  a,  b  = CoordArray ('a_ b_',   'NP-1', B-A           ).s() # [mm]
+  al, bl = ScalarArray('al bl',   'NP-1', a.magnitude() ).s() # [mm]
+  au, bu = CoordArray ('au bu',   'NP-1', a / al        ).s() # [1]
+
+  tan_theta = (au ^ bu) / (au & bu)     # [1]     bending
 
-  F = E + En * EFl
-  G = H + Hn * GHl
+  global mu, nu
+  mu, nu = CoordArray ('mu nu', 'NP-2', tan_theta ).s() # [1]
 
-  # diff_A_i = cost_ABCD.diff(vector_component(A, 0))
-  # global diff_A_i
-  # dbg('diff_A_i')
+  CostComponent('NP-3', sqnorm(mu - nu)) # [1]
+
+  dl2 = pow(al - bl, 2) # [mm^2]
+  CostComponent('NP-2', dl2 / (al*bl)) # [1]
+
+  # ---------- end of cost computation formulae ----------
 
   calculated = True
 
@@ -100,6 +169,12 @@ def cassign(v, assign_to, tmp_prefix):
   v = cse_prep_cprint(v, tmp_prefix)
   cprint(ourccode(v, assign_to=assign_to))
 
+def cassign_vector(v, assign_to, tmp_prefix):
+  ijk = 'i j k'.split(' ')
+  for ii in range(0, len(ijk)):
+    x = v & getattr(N, ijk[ii])
+    cassign(x, '%s[%d]' % (assign_to, ii), '%s_%s' % (tmp_prefix, ijk[ii]))
+
 def gen_diff(current, smalls):
   global j
   if not smalls:
@@ -145,46 +220,70 @@ def gen_point_index_macro(macro_basename, c_array_name, base_index):
 
 def gen_point_references():
   abcd = 'A B C D'.split(' ')
-  eh = 'E En H Hn'.split(' ')
-  for ehi in range(0, len(eh)):
-    gen_point_index_macro(eh[ehi], 'FIXED', ehi * 3)
 
+  gen_point_index_macro('E',  'INPUT', '3*0')
+  gen_point_index_macro('F0', 'INPUT', '3*1')
+  gen_point_index_macro('G0', 'INPUT', '3*(NP-2)')
+  gen_point_index_macro('H',  'INPUT', '3*(NP-1)')
+  cprintraw(         '#define NINPUT  ( 3*(NP-0) )')
+
+  gen_point_index_macro('En', 'PREP', '3*0')
+  gen_point_index_macro('Hn', 'PREP', '3*1')
+  cprintraw(         '#define NPREP   (3*2)')
+
+  cprintraw('#define NX_DIRECT 3*(NP-4)')
   cprint('#define POINT(PP) (')
-  cprint(' (PP) == 0   ? E :')
-  cprint(' (PP) == 1   ? F :')
-  cprint(' (PP) == N-2 ? G :')
-  cprint(' (PP) == N-1 ? H :')
-  cprint(' &X[((PP)-2)*3]')
+  cprint(' (PP) == 0    ? E :')
+  cprint(' (PP) == 1    ? F :')
+  cprint(' (PP) == NP-2 ? G :')
+  cprint(' (PP) == NP-1 ? H :')
+  cprint(' &X[3*((PP)-2)]')
   cprintraw(')')
 
+  cprintraw('#define EFlq X[ NX_DIRECT + 0 ]')
+  cprintraw('#define HGlq X[ NX_DIRECT + 1 ]')
+  cprintraw('#define NX    ( NX_DIRECT + 2 )')
+
   for ai in range(0, len(abcd)):
     cprintraw('#define %s POINT(P%+d)' % (abcd[ai], ai))
     gen_point_coords_macro(abcd[ai])
 
-def gen_calculate_FG():
-  cprint('#define CALCULAGE_F_G')
-  cassign(E,'E','tmp_E')
-  cassign(E,'F','tmp_F')
+  for si in iterations:
+    si.gen_references()
 
-def gen_calculate_cost():
-  cprint('#define CALCULATE_COST')
-  cassign(cost_ABCD,'P_cost','tmp_P_cost')
   cprintraw('')
 
-def gen_f_populate():
-  cprint('#define F_POPULATE')
-  cassign(cost_ABCD,'F','ftmp')
+def gen_prepare():
+  cprint('#define PREPARE')
+  cprint('memcpy(X, &INPUT[3*2], sizeof(double) * NX_DIRECT);')
+  for EH,EHs,FG0,FGs in ((E,'E', F0,'F'),
+                         (H,'H', G0,'G')):
+    EFHGv = FG0 - EH
+    EFHGl = EFHGv.magnitude()
+    EFHGlq = sqrt(EFHGl)
+    cassign_vector(EFHGv/EFHGl, EHs+'n', 'tmp_'+EHs)
+    cassign(EFHGlq, EHs+FGs+'lq', 'tmp_l'+EHs)
   cprintraw('')
 
-def gen_j_populate():
-  cprint('#define J_POPULATE')
-  gen_diff(result_dirnscaled, (sh*sh*la, th*th*la))
+def gen_calculate_FG():
+  cprintraw('#define DECLARE_F_G double F[3], G[3];')
+  cprint('#define CALCULATE_F_G')
+  cassign_vector(F,'F','tmp_F')
+  cassign_vector(G,'G','tmp_G')
+  cprintraw('')
+
+def gen_calculate_cost():
+  cprint('#define CALCULATE_COST')
+  cprint('double cost=0, P_cost;')
+  for si in iterations:
+    si.gen_calculate_cost()
   cprintraw('')
 
 def gen_C():
   gen_misc()
   gen_point_references()
-  #gen_calculate_FG()
+  gen_prepare()
+  gen_calculate_FG()
   gen_calculate_cost()
 
 def get_python():