X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?a=blobdiff_plain;f=symbolic.py;h=8680378ea06d906cbf58b4d230f09451e5ec281d;hb=c1001cc933c7a48935777908ba8338673cd091dc;hp=af5af8f1cbff32f3bd972f9f49976eab59761a16;hpb=51df2ada1c1d0da4ad56f4f5fdbe534b7008dfa0;p=moebius3.git diff --git a/symbolic.py b/symbolic.py index af5af8f..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 @@ -41,41 +47,102 @@ 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 * EFl - G = H + Hn * GHl + 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 - # diff_A_i = cost_ABCD.diff(vector_component(A, 0)) - # global diff_A_i - # dbg('diff_A_i') + 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 @@ -156,41 +223,50 @@ def gen_point_references(): gen_point_index_macro('E', 'INPUT', '3*0') gen_point_index_macro('F0', 'INPUT', '3*1') - gen_point_index_macro('G0', 'INPUT', '3*(N-2)') - gen_point_index_macro('H', 'INPUT', '3*(N-1)') - cprintraw( '#define NINPUT ( 3*(N-0) )') + 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 PREP_N (3*2)') + cprintraw( '#define NPREP (3*2)') - cprintraw('#define X_N_DIRECT 3*(N-4)') + 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(' (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 EFl X[ X_N_DIRECT + 0 ]') - cprintraw('#define GHl X[ X_N_DIRECT + 1 ]') - cprintraw('#define X_N ( X_N_DIRECT + 2 )') + 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) * X_N_DIRECT)') - cassign_vector((F0 - E).normalize(), 'En', 'tmp_En') - cassign_vector((G0 - H).normalize(), 'Hn', 'tmp_Hn') + 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') @@ -198,7 +274,9 @@ def gen_calculate_FG(): def gen_calculate_cost(): cprint('#define CALCULATE_COST') - cassign(cost_ABCD,'P_cost','tmp_P_cost') + cprint('double cost=0, P_cost;') + for si in iterations: + si.gen_calculate_cost() cprintraw('') def gen_C():