4
$\begingroup$

I've been bootstrapping myself to the very alien world of Mathematica and there came my first WTF moment:

RRP[P_, O_, Phi_, L_, M_] :=
  If[d <= L, O + l2*v0, {}] //. {
    v0 -> {Cos[Phi], Sin[Phi]},
    OP -> P - O,
    l1 -> Dot[OP, v0],
    v1 -> l1*v0,
    d  -> Norm[OP - v1],
    l2 -> l1 + M Sqrt[L^2 - d^2]
  };

What I wanted to achieve is to recursively apply these rules(which I felt similar to a context-free grammar) to form the desired expression. There was no cyclic dependency in the rules, nor were alternatives, so I kinda expected it to expand correctly. But it didn't!

Edit: I'm trying to calculate the position of the internal joint in a planar RRP linkage. Expressed in procedural language it's probably like this(pseudocode):

/* P: position of external R joint
 * O: reference point of P joint
 * Phi: direction of P joint
 * L: length of the R-R bar
 * M: +1/-1
 */
function RRP(P, O, Phi, L, M){
    OP=P-O;
    v0=[cos(Phi), sin(Phi)];
    l1=dot(OP, v0);
    v1=l1*v0;
    d=norm(OP-v1);
    l2=l1+M*sqrt(L^2-d^2);
    if(d<=L) return O+l2*v0;
    else return [];
}

Since I haven't figured out how to even debug here, I tried to output intermediate variables(or non-terminals in CFG terminology):

RRP[P_, O_, Phi_, L_, M_] :=
  {d, OP, v1, OP - v1} //. {
    v0 -> {Cos[Phi], Sin[Phi]},
    OP -> P - O,
    l1 -> Dot[OP, v0],
    v1 -> l1*v0,
    d -> Norm[OP - v1],
    l2 -> l1 + M Sqrt[L^2 - d^2]
    };
RRP[{0, 1}, {5, 0}, 0, Sqrt[2], 1]

And it gave me this:

{{6, Sqrt[26]}, {-5, 1}, {-5, 0}, {{0, -5}, {6, 1}}}

OP and v1 were as expected. But WTF was OP-v1? A Matrix?

I figure I can't do much with OP so I tried manually substituting v1 in the d-> rule:

RRP[P_, O_, Phi_, L_, M_] :=
  {d, OP, v1, OP - l1*v0} //. {
    v0 -> {Cos[Phi], Sin[Phi]},
    OP -> P - O,
    l1 -> Dot[OP, v0],
    v1 -> l1*v0,
    d -> Norm[OP - l1*v0],
    l2 -> l1 + M Sqrt[L^2 - d^2]
  };
RRP[{0, 1}, {5, 0}, 0, Sqrt[2], 1]

Now it gives me the correct answer:

{1, {-5, 1}, {-5, 0}, {0, 1}}

So now I'm VERY confused about this,

I don't even know if I'm following the proper/idiomatic way of constructing complex functions, thanks to the arcane Help that came with Mathematica... Could anyone please either point out the problem, or kindly enlighten me on how I could construct a equivalent function in "The Mathematica Way"?

EDIT: After some lookup I figured out that the evaluation order for OP - v1 was to blame. So I added Hold[] around the culprit and //ReleaseHold the whole thing after //. finishes building the expression. It works:

RRP[P_, O_, \[Phi]_, L_, M_] := (
    If[d <= L, O + l2*v0, {}] //. {
      v0 -> {Cos[\[Phi]], Sin[\[Phi]]},
      OP -> P - O,
      l1 -> Dot[OP, v0],
      v1 -> l1*v0,
      d -> Norm[Hold[OP - v1]],
      l2 -> l1 + M Sqrt[L^2 - d^2]
    }) // ReleaseHold;
RRP[{0, 1}, {5, 0}, 0, Sqrt[2], 1]

But I am still confused:

Could I know where I need Hold[] when I write the code, not as a workaround after I found out something wrong with it?

Improve this question
$\endgroup$
15
  • $\begingroup$ please delete the SO version. It is not a good idea to cross post $\endgroup$
    – rm -rf ♦
    Oct 4, 2012 at 14:52
  • $\begingroup$ Please also explain a bit of what you're trying to do here. $\endgroup$
    – Mr.Wizard
    Oct 4, 2012 at 14:52
  • $\begingroup$ @rm-rf okay I'm not quite familiar with the procedures here... I'll do right away $\endgroup$
    – smilekzs
    Oct 4, 2012 at 14:53
  • $\begingroup$ @Mr.Wizard should I write the expression in a procedural form(C-like)? $\endgroup$
    – smilekzs
    Oct 4, 2012 at 14:54
  • 2
    $\begingroup$ I'm guessing that you shouldn't use the single letter variable names C, D, E, I, K, N, or O. Best stick to lower-case for everything, as per Mr Wolfram $\endgroup$
    – cormullion
    Oct 4, 2012 at 15:35

2 Answers 2

Reset to default
14
$\begingroup$

While I can't follow your code, I guess your problem is caused by the fact that you get evaluation in between individual replacements, and the fact that Listable functions of several arguments (which includes operators like + and *) have quite peculiar behaviour. The fact that you get a matrix instead of a vector, as well as the fact that you can avoid it with a Hold point into this direction.

The problem with Listable functions and replacement

Consider the following two expressions

{1, 2} + 3
(*
==> {4, 5}
*)
{1, 2} + {3, 4}
(*
==> {4, 6}
*)

This is caused by the following rule:

  • If both arguments are lists, the operation is applied to pairs of corresponding elements.
  • If only one element is a list, then operation is applied to the elements of that list and the other argument.

Now consider the following expression:

{1, 2} + x
(*
==> {1 + x, 2 + x}
*)

Since x is not a list, the second bullet point above is used. This gives surprises if you replace x with a list:

{1, 2} + x /. x -> {3, 4}
(*
==> {{4, 5}, {5, 6}}
*)

This is because /. first evaluates the left hand side, and only then applies the rule. Therefore {1, 2} + x /. x -> {3, 4} first evaluates, according to the second bullet point, to {1 + x, 2 + x} /. x -> {3, 4} and only then does the replacement, leading to {1 + {3, 4}, 2 + {3, 4}} where for the individual sums again the second bullet point rule applies, to give {{1+3, 1+4}, {2+3, 2+4}} which gives the result above. Note how in the process, you've generated a matrix from your vectors.

Interaction with //.

Now consider the following two expressions:

x + y //. {x -> {1,2}, y -> {3,4}}
(*
==> {4, 6}
*)
x + y //. {x -> {1,2}, y -> z, z -> {3,4}}
(*
==> {{4, 5}, {5, 6}}
*)

As you see, the first one gives the expected expression, while the second one doesn't, despite the fact that looking at the replacements, it apparently should be equivalent.

The reason for the difference is that after each round of replacement, the expression is evaluated again, before the next round of replacements is done. So in the first case, after replacing x and y, you end up with {1, 2} + {3, 4} which evaluates to {4, 6} whithout any further replacements. However in the second case, the first round of replacement leaves you with {1, 2} + z which, as described above, will evaluate to {1 + z, 2 + z}, and that is what the second round of replacements acts on, replacing z with {3, 4}.

Solution

If you just don't want to evaluate anything before all replacement rules have been applied, just wrap the original expression in Hold and the whole expression including the replacement rules in ReleaseHold:

ReleaseHold[Hold[x + y] //. {x -> {1,2}, y -> z, z -> {3,4}}]
(*
==> {4, 6}
*)

If you want the original expression to be evaluated just once before applying the rules, one trick is to put it in an extra initial rule which is not deferred (i.e. ->, not :>):

f[z] = x;
ReleaseHold[Hold[initial] //.
  {initial -> f[z] + y, x -> {1,2}, y -> z, z -> {3,4}}]
(*
==> {4, 6}
*)

Without this initial rule trick you'd get

ReleaseHold[Hold[f[z] + y] //. {x -> {1,2}, y -> z, z -> {3,4}}]
(*
==> {3 + f[{3, 4}], 4 + f[{3, 4}]}
*)

because z would be replaced before f got applied on it.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Yes exactly(although I figured out the vector-scalar relationship earlier myself earlier). You made it clear without doubt. $\endgroup$
    – smilekzs
    Oct 4, 2012 at 16:07
2
$\begingroup$

After staring at this for a bit I think I understand what you're trying to do, but I'm not sure what I have to say will be satisfactory. It seems that you are giving a series of transformation rules that you want to be carried out in isolation, that is apart from the normal evaluation process. For this reason Hold is naturally needed, and it should wrap the entire expression, not just one element.

RRP[P_, O_, Phi_, L_, M_] :=
 ReleaseHold[
  Hold[If[d <= L, O + l2*v0, {}]] //.
   {v0 -> {Cos[Phi], Sin[Phi]},
    OP -> P - O,
    l1 -> Dot[OP, v0],
    v1 -> l1*v0,
    d -> Norm[OP - v1],
    l2 -> l1 + M Sqrt[L^2 - d^2]}
 ]

If this technique does not work in another case of the same method please post that case as well.

Of course it's probably easier in this case just to translate your pseudocode literally: 

RRP[P_, O_, Phi_, L_, M_] :=
 Module[{OP, l1, l2, v0, v1, d},
  OP = P - O;
  v0 = {Cos[Phi], Sin[Phi]};
  l1 = OP.v0;
  v1 = l1*v0;
  d = Norm[OP - v1];
  l2 = l1 + M*Sqrt[L^2 - d^2];
  If[d <= L, O + l2*v0, {}]
 ]

RRP[{0, 1}, {5, 0}, 0, Sqrt[2], 1]

{1, 0}

I don't believe I properly understand your program yet, but I'll show you how I might go about visualizing what is taking place.

We can use FixedPointList and /. in place of //., which will allow us to keep the steps of evaluation. Using that we can compare the versions with and without Hold:

ClearAll[RRP]

RRP[P_, O_, \[Phi]_, L_, M_] := 
  FixedPointList[# /. {v0 -> {Cos[\[Phi]], Sin[\[Phi]]}, OP -> P - O, 
       l1 -> Dot[OP, v0], v1 -> l1*v0, d -> Norm[Hold[OP - v1]], 
       l2 -> l1 + M Sqrt[L^2 - d^2]} &, If[d <= L, O + l2*v0, {}]] // 
   ReleaseHold

RRP[{0, 1}, {5, 0}, 0, Sqrt[2], 1] // Column

Output:

If[d<=Sqrt[2],{5,0}+l2 v0,{}]

If[Norm[OP-v1]<=Sqrt[2],{5,0}+(Sqrt[2-d^2]+l1) {1,0},{}]

If[Sqrt[Abs[-5-l1 v0]^2+Abs[1-l1 v0]^2]<=Sqrt[2],{5,0}+(Sqrt[2-Norm[OP-v1]^2]+OP.v0) {1,0},{}]

If[Sqrt[1+Abs[-5-OP.v0]^2]<=Sqrt[2],{5,0}+(Sqrt[2-Norm[{-5,1}-l1 v0]^2]+{-5,1}.{1,0}) {1,0},{}]

{Sqrt[1-Abs[-5-OP.v0]^2],0}

{1,0}

{1,0}

Now without Hold:

ClearAll[RRP]

RRP[P_, O_, \[Phi]_, L_, M_] := 
 FixedPointList[# /. {v0 -> {Cos[\[Phi]], Sin[\[Phi]]}, OP -> P - O, 
     l1 -> Dot[OP, v0], v1 -> l1*v0, d -> Norm[OP - v1], 
     l2 -> l1 + M Sqrt[L^2 - d^2]} &, If[d <= L, O + l2*v0, {}]]

RRP[{0, 1}, {5, 0}, 0, Sqrt[2], 1] // Column

Output:

If[d<=Sqrt[2],{5,0}+l2 v0,{}]

If[Norm[OP-v1]<=Sqrt[2],{5,0}+(Sqrt[2-d^2]+l1) {1,0},{}]

If[Sqrt[Abs[-5-l1 v0]^2+Abs[1-l1 v0]^2]<=Sqrt[2],{5,0}+(Sqrt[2-Norm[OP-v1]^2]+OP.v0) {1,0},{}]

If[{Sqrt[Abs[-5-OP.v0]^2+Abs[1-OP.v0]^2],Sqrt[26]}<=Sqrt[2],{5,0}+(Sqrt[2-Norm[{-5,1}-l1 v0]^2]+{-5,1}.{1,0}) {1,0},{}]

If[{6,Sqrt[26]}<=Sqrt[2],{5,0}+(Sqrt[2-Norm[{-5,1}-OP.v0 {1,0}]^2]+{-5,1}.{1,0}) {1,0},{}]

If[{6,Sqrt[26]}<=Sqrt[2],{5,0}+(Sqrt[2-Norm[{-5,1}-{-5,1}.{1,0} {1,0}]^2]+{-5,1}.{1,0}) {1,0},{}]

If[{6,Sqrt[26]}<=Sqrt[2],{5,0}+(Sqrt[2-Norm[{-5,1}-{-5,1}.{1,0} {1,0}]^2]+{-5,1}.{1,0}) {1,0},{}]

Perhaps you can see from this where things go wrong.

Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ As I told the OP, using the protected symbol O as an argument might not be the best idea... $\endgroup$
    – J. M.'s persistent exhaustion ♦
    Oct 4, 2012 at 15:13
  • 1
    $\begingroup$ @J.M. I am not suggesting it! I am trying to show a method of debugging; I don't what to change the code more than necessary. $\endgroup$
    – Mr.Wizard
    Oct 4, 2012 at 15:14
  • $\begingroup$ Yes I have already figured out that OP-v1 expanded first to {_,_}-v1 then {_-v1,_-v1}, before v1 even got expanded, where the - was understood as operation between vector and scalar, all while I really mean to "expand both OP and v1 before applying the vector difference"... $\endgroup$
    – smilekzs
    Oct 4, 2012 at 15:15
  • $\begingroup$ Hmm I think I get the idea about how and why Hold[] works now. I don't really like the idea of injecting syntactic noise into mathematical notation, so could you please suggest a way I could rewrite the code in a more idiomatic way? $\endgroup$
    – smilekzs
    Oct 4, 2012 at 15:29
  • 1
    $\begingroup$ Oh... So I could just write procedural code as I'm used to. I thought it has to be purely functional! $\endgroup$
    – smilekzs
    Oct 4, 2012 at 15:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.