title: Recursively calling functions within functions in Python (trying to replicate MATLAB behaviour)
body: <p>In MATLAB this function (by Hao Zhang) calls itself</p>
  
  <pre><code>function r=rotmat2expmap(R)
  
  % Software provided by Hao Zhang
  % http://www.cs.berkeley.edu/~nhz/software/rotations
  
  r=quat2expmap(rotmat2quat(R));
  </code></pre>
  
  <p>as an argument to the function </p>
  
  <pre><code>function [r]=quat2expmap(q)
  
  % Software provided by Hao Zhang
  % http://www.cs.berkeley.edu/~nhz/software/rotations
  %
  % function [r]=quat2expmap(q)
  % convert quaternion q into exponential map r
  % 
  % denote the axis of rotation by unit vector r0, the angle by theta
  % q is of the form (cos(theta/2), r0*sin(theta/2))
  % r is of the form r0*theta
  
    if (abs(norm(q)-1)&gt;1E-3)
      error('quat2expmap: input quaternion is not norm 1');
    end
    sinhalftheta=norm(q(2:4));
    coshalftheta=q(1);
    r0=q(2:4)/(norm(q(2:4))+eps);
    theta=2*atan2(sinhalftheta,coshalftheta);
    theta=mod(theta+2*pi,2*pi);
    %if (theta&gt;pi), theta=2*pi-theta; r0=-r0; end
    if (theta&gt;pi)
      theta=2*pi-theta; 
      r0=-r0; 
    end
    r=r0*theta;
  </code></pre>
  
  <p>Now if we pass a rotation matrix to the first function something along the lines of</p>
  
  <pre><code>R =
  
      0.9940    0.0773   -0.0773
     -0.0713    0.9945    0.0769
      0.0828   -0.0709    0.9940
  </code></pre>
  
  <p>It recursively calculates the correct result (in this case) simply:</p>
  
  <pre><code>r =
  
     -0.0741   -0.0803   -0.0745
  </code></pre>
  
  <p>Alas this is in MATLAB and it works fine (the original author knew what he was doing). I have not quite managed to get the same functionality to work in Python (I am effectively translating the code), I am going wrong somewhere:</p>
  
  <pre><code>def rotmat2expmap(R):
      """
      Convert rotation matrix to exponential mapping.
      Based on G.W. Taylor's MATLAB equivalent.
      """
      r = quat2expmap(rotmat2expmap(R))
      return r
  
  def quat2expmap(q):
      """Converts quaternion q (rotation matrix) into exponential map r.
      Provided by Hao Zhang and G.W. Taylor.
  
      Denote the axis of rotation by unit vector r0, the angle by theta
      q is of the form (cos(theta/2), r0*sin(theta/2))
      r is of the form r0*theta
      """
      if abs(np.linalg.norm(q,2)-1) &gt; 1e-3:
          print('quat2expmap: input quaternion is not norm 1')
  
      # Implement to simulate MATLAB like linear array structure
      temp = q.T.flatten()
  
      sinhalftheta = np.linalg.norm(temp[1:4],2)
      coshalftheta = temp[0]
      r0 = temp[1:4]/(np.linalg.norm(temp[1:4],2) + np.spacing(1))
      theta = 2*math.atan2(sinhalftheta,coshalftheta)
      theta = fmod(theta+2*pi,2*pi) # Remainder after division (modulo operation)
      if theta &gt; pi:
          theta = 2*pi-theta
          r0 = -r0
      r = r0*theta
  
      return r
  </code></pre>
  
  <p>If I try to run this (with the same example R) then the number of loops maxes out and the whole thing crashes. </p>
  
  <p>Anyone got any fancy ideas?</p>