So I think it has become clear what I'm looking to do. There are many levels of subexpressions, and equations may be defined by many different subexpressions. For those of you familiar with classical mechanics, you are probably aware of the explosion of the number of terms that arises as the number of rigid bodies increases. The following is a snippet of code from an Autolev script that I've used to formulate the equations of motion for a benchmark Whipple bicycle model, with the interpreter response being indicated by an '->' and the introduction of intermediate variables beginning with the letter 'Z':
(98) %--------------------------------------------------------------------- % (99) % ANGULAR RELATIONSHIPS % (100) %--------------------------------------------------------------------- % (101) % FRAME YAW (102) SIMPROT(N,YF,3,q3) -> (103) Z1 = COS(q3) -> (104) Z2 = SIN(q3) -> (105) N_YF = [Z1, -Z2, 0; Z2, Z1, 0; 0, 0, 1] (106) % FRAME LEAN (107) SIMPROT(YF,LF,1,q4) -> (108) Z3 = COS(q4) -> (109) Z4 = SIN(q4) -> (110) YF_LF = [1, 0, 0; 0, Z3, -Z4; 0, Z4, Z3] (111) % REAR WHEEL ROTATION (112) SIMPROT(LF,R,2,q5) -> (113) Z5 = COS(q5) -> (114) Z6 = SIN(q5) -> (115) LF_R = [Z5, 0, Z6; 0, 1, 0; -Z6, 0, Z5] (116) % FRAME PITCH (117) SIMPROT(LF,B,2,q6) -> (118) Z7 = COS(q6) -> (119) Z8 = SIN(q6) -> (120) LF_B = [Z7, 0, Z8; 0, 1, 0; -Z8, 0, Z7] (121) % INTERMEDIATE STEER FRAME (122) SIMPROT(B,SF,2,lambda) -> (123) Z9 = COS(lambda) -> (124) Z10 = SIN(lambda) -> (125) B_SF = [Z9, 0, Z10; 0, 1, 0; -Z10, 0, Z9] (126) % STEERING ANGLE (127) SIMPROT(SF,H,3,q7) -> (128) Z11 = COS(q7) -> (129) Z12 = SIN(q7) -> (130) SF_H = [Z11, -Z12, 0; Z12, Z11, 0; 0, 0, 1] (131) % HANDLEBAR MASS FRAME ROTATION (132) SIMPROT(H,HF,2,-lambda) -> (133) H_HF = [Z9, 0, -Z10; 0, 1, 0; Z10, 0, Z9] (134) % FRONT WHEEL ROTATION (135) SIMPROT(H,F,2,q8) -> (136) Z13 = COS(q8) -> (137) Z14 = SIN(q8) -> (138) H_F = [Z13, 0, Z14; 0, 1, 0; -Z14, 0, Z13] (139) %--------------------------------------------------------------------- % (140) % POSITION VECTORS (141) %--------------------------------------------------------------------- % (142) % Locate the rear wheel center of mass (143) P_NO_RO> = q1*N1> + q2*N2> - rrt*N3> - rr*LF3> -> (144) P_NO_RO> = -rr*LF3> + q1*N1> + q2*N2> - rrt*N3> (145) % Rear contact points (146) % We need two points here because one is fixed to the wheel, the (147) % other is fixed to the ground. Later when we need to impose the (148) % nonholonomic constraint, we form the velocity by treating RN as (149) % a point FIXED in C (150) P_RO_NR> = rr*LF3> + rrt*N3> -> (151) P_RO_NR> = rr*LF3> + rrt*N3> (152) P_RO_RN> = P_RO_NR> -> (153) P_RO_RN> = rr*LF3> + rrt*N3> (154) % Locate the frame center of mass (155) P_RO_BO> = d4*B1> + d5*B3> -> (156) P_RO_BO> = d4*B1> + d5*B3> (157) % Locate the front wheel center of mass (158) P_RO_FO> = d1*SF1> + d3*SF3> + d2*H1> -> (159) P_RO_FO> = d2*H1> + d1*SF1> + d3*SF3> (160) % Locate the fork center of mass (161) P_RO_HO> = d1*SF1> + d6*HF1> + d7*HF3> -> (162) P_RO_HO> = d6*HF1> + d7*HF3> + d1*SF1> (163) % Front contact points (164) % We need two points here because one is fixed to the wheel, the (165) % other is fixed to the ground. Later when we need to impose the (166) % nonholonomic constraint, we form the velocity by treating FN as (167) % a point FIXED in F (168) P_FO_NF> = rf*UNITVEC(N3>-DOT(H2>,N3>)*H2>) + rft*N3> -> (169) Z15 = Z9*Z11 -> (170) Z16 = Z10*Z11 -> (171) Z17 = Z9*Z12 -> (172) Z18 = Z10*Z12 -> (173) H_B = [Z15, Z12, -Z16; -Z17, Z11, Z18; Z10, 0, Z9] -> (174) Z19 = Z4*Z8 -> (175) Z20 = Z3*Z8 -> (176) Z21 = Z4*Z7 -> (177) Z22 = Z3*Z7 -> (178) B_YF = [Z7, Z19, -Z20; 0, Z3, Z4; Z8, -Z21, Z22] -> (179) Z23 = Z1*Z7 - Z2*Z19 -> (180) Z24 = Z1*Z19 + Z2*Z7 -> (181) Z25 = Z2*Z3 -> (182) Z26 = Z1*Z3 -> (183) Z27 = Z1*Z8 + Z2*Z21 -> (184) Z28 = Z2*Z8 - Z1*Z21 -> (185) B_N = [Z23, Z24, -Z20; -Z25, Z26, Z4; Z27, Z28, Z22] -> (186) Z29 = Z15*Z23 - Z12*Z25 - Z16*Z27 -> (187) Z30 = Z12*Z26 + Z15*Z24 - Z16*Z28 -> (188) Z31 = Z4*Z12 - Z15*Z20 - Z16*Z22 -> (189) Z32 = Z18*Z27 - Z11*Z25 - Z17*Z23 -> (190) Z33 = Z11*Z26 + Z18*Z28 - Z17*Z24 -> (191) Z34 = Z4*Z11 + Z17*Z20 + Z18*Z22 -> (192) Z35 = Z9*Z27 + Z10*Z23 -> (193) Z36 = Z9*Z28 + Z10*Z24 -> (194) Z37 = Z9*Z22 - Z10*Z20 -> (195) H_N = [Z29, Z30, Z31; Z32, Z33, Z34; Z35, Z36, Z37] -> (196) Z38 = 1 - Z34^2 -> (197) Z39 = Z38^0.5 -> (198) Z40 = 1/Z39 -> (199) Z41 = Z34/Z39 -> (200) Z42 = rf*Z41 -> (201) Z43 = rft + rf*Z40 -> (202) P_FO_NF> = -Z42*H2> + Z43*N3> Comments: Lines 98-138: Here is where the angular relationship of the various rigid bodies and intermediate frames are defined. The 'SIMPROT' command stands for Simple Rotation, and it essentially just generate a rotation matrix, but notice that it automatically introduces intermediate variables for the sin and cosine of the angle of rotation. Lines 139-162: The notation here is a little funky, but basically 'P_NO_RO>' denotes the position vector (hence the 'P'), from the point NO, to the point RO. Autolev recognizes 'NO' as the origin of a Newtonian reference frame that I declared earlier (not shown), and 'RO' as the center of mass location of the rigid body 'R' (also declared earlier). Anything with a '>' at the end of it denotes that it is a vector, and most of the vectors here are the body fixed vectors of each rigid body (they are created automatically upon declaration of a rigid body). Once the orientations of all of the bodies has been declared, you can work with whichever coordinate system is most convenient when describing the location of a center of mass or a point. Lines 163-202: This best illustrates the behavior I am trying to achieve. Upon entering the definition of 'P_FO_NF>' (the position vector from the center of the front wheel to the front wheel - ground contact point), the Autolev interpreter generates 28 intermediate 'Z' variables before finally forming the final expression. All of these intermediate variables somehow involve the generalized coordinates of the problem at hand (most happen to be angles) and hence they would vary throughout the course of a numerical integration. For comparison, if P_FO_NF> is expressed in the Newtonian frame and no substitution is used, it looks like: -> (224) P_FO_NF> = rf*(SIN(q4)*COS(q7)+SIN(q7)*COS(q4)*SIN(lambda +q6))*(SIN(q7) *COS(q3)*COS(lambda+q6)+SIN(q3)*(COS(q4)*COS(q7)- SIN(q4)*SIN(q7)*SIN(lambda+q6)) )/(1-(SIN(q4)*COS(q7)+SIN(q7)*COS(q4)*SIN(lambda+q6))^2)^0.5*N1> + rf*(SIN(q4)*C OS(q7)+SIN(q7)*COS(q4)*SIN(lambda+q6))*(SIN(q3)*SIN(q7)*COS(lambda+q6)- COS(q3)*( COS(q4)*COS(q7)-SIN(q4)*SIN(q7)*SIN(lambda+q6)))/(1- (SIN(q4)*COS(q7)+SIN(q7)*COS (q4)*SIN(lambda+q6))^2)^0.5*N2> + (rft+rf*(1- (SIN(q4)*COS(q7)+SIN(q7)*COS(q4)*SI N(lambda+q6))^2)^0.5)*N3> It is easy to see that there are many repeated quantities in the above expression and it makes sense that they only be computed once. When forming expressions for the velocities and accelerations, one can imagine how long the expressions get after the application of the chain rule.... many pages long. I don't know if the method Autolev is identifying and creating intermediate variables is the best way to do it because I have seen cases that don't seem to actually simplify the equations much, but for the most part, it seems to work well. And, as was mentioned by somebody else, there must be a trade off between memory usage and flop usage. I the problems I have worked on, using extra memory used has been worth the savings in flops by a large factor. For example, the equations of motion for the Whipple bicycle generated by Autolev are about 77kb, while those same equations generated by another author's method are about 3.5mb, almost 2 orders of magnitude larger. Needless to say, numerical integration times are much faster with the smaller more compact representation of the equations of motion. Essentially what I need to do is parse every equation and identify subexressions (and then parse those subexpression and identify subsubexpressions....), and introduce variables for any that are used more than once. If a subexpression is found, but doesn't occur more than once in any of the preceding equations, but then later shows up in another equation, it would then be introduced as an intermediate variable. If they aren't used more than once, then it doesn't make sense to introduce an intermediate variable for them, unless the goal is simply just readability. Any ideas on how to start parsing each equation and accumulating a list of subexpressions? I haven't had time to start coding this up in python yet, hopefully I'll get started on this over the summer. Thanks, ~Luke On Jun 15, 11:28 pm, "Robert Kern" <[EMAIL PROTECTED]> wrote: > On Sun, Jun 15, 2008 at 17:46, Ondrej Certik <[EMAIL PROTECTED]> wrote: > > Anyway, I created an issue for this: > > >http://code.google.com/p/sympy/issues/detail?id=891 > > I've attached some nominally working code to the issue. Playing with > it, I see some suboptimal results. First, it's quadratic, but I'm not > clever enough to think of a way around that. Second, sympy's default > representation of expressions isn't well-suited for this algorithm (or > the other way around). For example: > > In [1]: from cse import cse > > In [2]: cse((x-y)*(z-y) + sqrt((x-y)*(z-y))) > Out[2]: > ⎛ ⎽⎽⎽⎽ ⎽⎽⎽⎽⎞ > ⎝[(x₀, -y), (x₁, x + x₀), (x₂, x₀ + z)], x₁*x₂ + ╲╱ x₁ *╲╱ x₂ ⎠ > > We see that -y gets picked up as term that should be pulled out > although it only shows up in the context of (x-y). For the typical use > case, subtraction is an atomic operation, and the representation of > (x-y) as Add(Mul(NegativeOne(-1), Symbol('y')), Symbol('x')) gets in > the way. > > Although the input expression has both x1 and x2 under the same > sqrt(), it gets broken out before the cse() function gets to look at > it. It would be nice to stuff everything in that term under the same > sqrt() before cse() operates. This kind of thing pops up from time to > time in my experience. > > So I think there needs to be a bit of a preprocessing step to optimize > the expression specifically for cse() before it does its magic. > > Generalizing this to a list of expressions (e.g. my old use case) is > an exercise left to the reader. > > -- > Robert Kern > > "I have come to believe that the whole world is an enigma, a harmless > enigma that is made terrible by our own mad attempt to interpret it as > though it had an underlying truth." > -- Umberto Eco --~--~---------~--~----~------------~-------~--~----~ You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to sympy@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sympy?hl=en -~----------~----~----~----~------~----~------~--~---
