X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?a=blobdiff_plain;f=symbolic.py;h=8680378ea06d906cbf58b4d230f09451e5ec281d;hb=c1001cc933c7a48935777908ba8338673cd091dc;hp=e44846eda32231448a75212d992cde586fdeb8ae;hpb=553d1b0eb2931fb104c9cdaa79c7c391595485a6;p=moebius3.git diff --git a/symbolic.py b/symbolic.py old mode 100755 new mode 100644 index e44846e..8680378 --- a/symbolic.py +++ b/symbolic.py @@ -1,188 +1,297 @@ -#!/usr/bin/python3 +from __future__ import print_function + +from sympy.vector.vector import * from sympy import * import itertools +from sympy.utilities.lambdify import lambdify, implemented_function +import numpy as np -import sys - -import sys, codecs -if sys.stdout.encoding is None: - sys.stdout = codecs.open("/dev/stdout", "w", 'utf-8') - -init_printing(use_unicode=True, num_columns=280) +from moedebug import * -r, theta, s, la, mu, kappa = symbols('r theta s lambda mu kappa') +# When cse() is called for something containing BaseVector, it +# produces infinite recursion. +#def cse(x, *a, **kw): return ((), (x,)) -# start original formulation -# rightvars replaces - -p_start = Matrix([ - r * (1 - cos(theta)), - r * sin(theta), - mu * s, -]) - -p_rightvars = p_start.subs( theta, s/r ).subs( r, 1/la ) +def dprint(*args): + if not dbg_enabled(): return + print(*args) def dbg(*args): + if not dbg_enabled(): return for vn in args: print('\n ' + vn + '\n') pprint(eval(vn)) print('\n =\n') pprint(cse(eval(vn))) -dbg('p_rightvars') - -p_dirn_rightvars = diff(p_rightvars, s) - -dbg('p_dirn_rightvars') - -zeta = Wild('zeta') - -p_nosing = (p_rightvars - .replace( 1-cos(zeta) , 2*sin(zeta/2)**2 ) - .replace( sin(zeta)**2 , zeta*sinc(zeta)*sin(zeta) ) - ) -p_nosing[1] = (p_nosing[1] - .replace( sin(zeta) , zeta * sinc(zeta) ) - ) - -dbg('p_nosing') - -t = symbols('t') - -q_owncoords = p_nosing.replace(s,t).replace(la,-la) -q_dirn_owncoords = p_dirn_rightvars.replace(s,t).replace(la,-la) - -dbg('q_owncoords','q_dirn_owncoords') -dbg('q_owncoords.replace(t,0)','q_dirn_owncoords.replace(t,0)') - -p2q_translate = p_nosing -#p2q_rotate_2d = Matrix([ p_dirn_rightvars[0:2], - -#p2q_rotate = eye(3) -#p2q_rotate[0:2, 0] = Matrix([ p_dirn_rightvars[1], -p_dirn_rightvars[0] ]) -#p2q_rotate[0:2, 1] = p_dirn_rightvars[0:2] - -p2q_rotate = Matrix([[ cos(theta), sin(theta), 0 ], - [ -sin(theta), cos(theta), 0 ], - [ 0 , 0, 1 ]]).subs(theta,la*s) -#p2q_rotate.add_col([0,0]) -#p2q_rotate.add_row([0,0,1]) - -dbg('p2q_rotate') - -q_dirn_maincoords = p2q_rotate * q_dirn_owncoords; -q_maincoords = p2q_rotate * q_owncoords + p2q_translate - -dbg('diff(p_dirn_rightvars,s)') -dbg('diff(q_dirn_maincoords,t)') -dbg('diff(q_dirn_maincoords,t).replace(t,0)') - -assert(Eq(p2q_rotate * Matrix([0,1,mu]), p_dirn_rightvars)) - -#for v in 's','t','la','mu': -# dbg('diff(q_maincoords,%s)' % v) - -#print('\n eye3 subs etc.\n') -#dbg('''Eq(eye(3) * Matrix([1,0,mu]), -# p_dirn_rightvars .cross(Matrix([0,0,1]) .subs(s,0)))''') - -#dbg('''Eq(p2q_rotate * Matrix([1,0,mu]), -# p_dirn_rightvars .cross(Matrix([0,0,1])))''') - -#eq = Eq(qmat * q_dirn_owncoords_0, p_dirn_rightvars) -#print -#pprint(eq) -#solve(eq, Q) - -dbg('q_maincoords.replace(t,0)','q_dirn_maincoords.replace(t,0)') - -dbg('q_maincoords','q_dirn_maincoords') - -sinof_mu = sin(atan(mu)) -cosof_mu = cos(atan(mu)) - -dbg('cosof_mu','sinof_mu') - -o2p_rotate1 = Matrix([[ 1, 0, 0 ], - [ 0, cosof_mu, +sinof_mu ], - [ 0, -sinof_mu, cosof_mu ]]) - -check_dirn_p_s0 = o2p_rotate1 * p_dirn_rightvars.replace(s,0) -check_dirn_p_s0.simplify() -dbg('check_dirn_p_s0') - -o2p_rotate2 = Matrix([[ cos(kappa), 0, -sin(kappa) ], - [ 0, 1, 0 ], - [ +sin(kappa), 0, cos(kappa) ]]) - -p_dirn_orgcoords = o2p_rotate2 * o2p_rotate1 * p_dirn_rightvars - -check_dirn_p_s0 = p_dirn_orgcoords.replace(s,0) -check_dirn_p_s0.simplify() -dbg('check_dirn_p_s0') - -check_accel_p_s0 = diff(p_dirn_orgcoords,s).replace(s,0) -check_accel_p_s0.simplify() -dbg('check_accel_p_s0') - -q_dirn_orgcoords = o2p_rotate2 * o2p_rotate1 * q_dirn_maincoords; -q_orgcoords = o2p_rotate2 * o2p_rotate1 * q_maincoords; -dbg('q_orgcoords','q_dirn_orgcoords') - -sh, th = symbols('alpha beta') - -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) - -print('----------------------------------------') -dbg('q_sqparm', 'q_dirn_sqparm') -print('----------------------------------------') -for v in 'sh','th','la','mu': - dbg('diff(q_sqparm,%s)' % v) - dbg('diff(q_dirn_sqparm,%s)' % v) -print('----------------------------------------') - -gamma = symbols('gamma') - -q_dirn_dirnscaled = q_dirn_sqparm * gamma - -result_dirnscaled = q_sqparm.col_join(q_dirn_dirnscaled) -dbg('result_dirnscaled') - -params = ('sh','th','la','mu','gamma','kappa') +def sqnorm(v): return v & v + +N = CoordSysCartesian('N') + +calculated = False + +def vector_symbols(vnames): + out = [] + for vname in vnames.split(' '): + v = Vector.zero + for cname in 'i j k'.split(' '): + v += getattr(N, cname) * symbols(vname + '_' + cname) + out.append(v) + return out + +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') + +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 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 + + # ---------- 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 + + global mu, nu + mu, nu = CoordArray ('mu nu', 'NP-2', tan_theta ).s() # [1] + + 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 + +def ourccode(*a, **kw): + return ccode(*a, user_functions={'sinc':'sinc'}, **kw) + +def cprintraw(*s): + print(*s) + +def cprint(s): + for l in s.split('\n'): + cprintraw(l, '\\') + +def cse_prep_cprint(v, tmp_prefix): + # => v, but also having cprint'd the common subexpression assignments + sym_iter = map((lambda i: symbols('%s%d' % (tmp_prefix,i))), + itertools.count()) + (defs, vs) = cse(v, symbols=sym_iter) + for defname, defval in defs: + cprint('double '+ourccode(defval, assign_to=defname)) + return vs[0] + +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: + j = zeros(len(params),0) for param in params: global d paramv = eval(param) d = diff(current, paramv) dbg('d') - print(ccode(d, assign_to='FOO')) - (defs, v) = cse(d, symbols=map((lambda i: symbols('tmp%d' % i)), - itertools.count())) - for vn, val in defs: - print(ccode(val, assign_to=vn)) - print(ccode(v[0], assign_to='BAR')) + j = j.row_join(d) + dbg('j') + j = cse_prep_cprint(j, 'jtmp') + for ix in range(0, j.cols): + cprint(ourccode(j.col(ix), 'J_COL')) + cprint('J_END_COL(%d)' % ix) else: small = smalls[0] smalls = smalls[1:] - print('if ',small,' nonsmall {') + cprint('if (!IS_SMALL(' + ourccode(small) + ')) {') gen_diff(current, smalls) - print('} else { /* %s small */' % small) + cprint('} else { /* %s small */' % small) gen_diff(current.replace( sinc(small), 1 - small*small/factorial(3) - small**4/factorial(5), ), smalls ) - print('} /* %s small */' % small) - -gen_diff(result_dirnscaled, (sh*sh*la, th*th*la)) - -#bad = q_orgcoords[0] -#badd = diff(bad, la) - -#dbg('bad','badd') + cprint('} /* %s small */' % small) + +def gen_misc(): + cprintraw('// AUTOGENERATED - DO NOT EDIT\n') + +def gen_point_coords_macro(macro_basename): + ijk = 'i j k'.split(' ') + for ii in range(0, len(ijk)): + cprintraw('#define %s_%s (%s[%d])' + % (macro_basename, ijk[ii], macro_basename, ii)) + +def gen_point_index_macro(macro_basename, c_array_name, base_index): + cprintraw('#define %s (&%s[%s])' + % (macro_basename, c_array_name, base_index)) + gen_point_coords_macro(macro_basename) + +def gen_point_references(): + abcd = 'A B C D'.split(' ') + + 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) == 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]) + + for si in iterations: + si.gen_references() + + cprintraw('') + +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_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_prepare() + gen_calculate_FG() + gen_calculate_cost() + +def get_python(): + # https://github.com/sympy/sympy/issues/13642 + # "lambdify sinc gives wrong answer!" + out = q_sqparm + sinc_fixed = Function('sinc_fixed') + implemented_function(sinc_fixed, lambda x: np.sinc(x/np.pi)) + out = out.subs(sinc,sinc_fixed) + p = list(map(eval,params)) + return lambdify(p, out)