Source Code for Module pylal.sphericalutils

  1  # Copyright (C) 2010  Nickolas Fotopoulos 
  2  # 
  3  # This program is free software; you can redistribute it and/or modify it 
  4  # under the terms of the GNU General Public License as published by the 
  5  # Free Software Foundation; either version 2 of the License, or (at your 
  6  # option) any later version. 
  7  # 
  8  # This program is distributed in the hope that it will be useful, but 
  9  # WITHOUT ANY WARRANTY; without even the implied warranty of 
 10  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General 
 11  # Public License for more details. 
 12  # 
 13  # You should have received a copy of the GNU General Public License along 
 14  # with this program; if not, write to the Free Software Foundation, Inc., 
 15  # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA. 
 16  """ 
 17  Utilities to perform operations on (polar, azimuthal) vectors. 
 18  """ 
 19   
 20  from __future__ import division 
 21   
 22  import numpy as np 
 23  np.seterr(all="raise") 
 24   
 25  import lal 
 26   
 27  __author__ = "Nickolas Fotopoulos <nvf@gravity.phys.uwm.edu>" 
 28   
 29   
 30  # 
 31  # Rotation utilities 
 32  # 
 33   
 34 -def rotate_euler(sph_coords, alpha, beta, gamma): 
 35      """ 
 36      Take an Nx2 array of N (theta, phi) vectors on the unit sphere 
 37      (that is, (polar, azimuthal) angles in radians) and apply 
 38      rotations through the Euler angles alpha, beta, and gamma 
 39      in radians, using the ZXZ convention. 
 40      """ 
 41      c = np.cos 
 42      s = np.sin 
 43   
 44      # Define rotation matrix 
 45      R = np.array( 
 46          [[c(alpha) * c(gamma) - c(beta) * s(alpha) * s(gamma), 
 47            c(gamma) * s(alpha) + c(alpha) * c(beta) * s(gamma), 
 48            s(beta) * s(gamma)], 
 49           [-c(beta) * c(gamma) * s(alpha) - c(alpha) * s(gamma), 
 50            c(alpha) * c(beta) * c(gamma) - s(alpha) * s(gamma), 
 51            c(gamma) * s(beta)], 
 52           [s(alpha) * s(beta), -c(alpha) * s(beta), c(beta)]], dtype=float) 
 53   
 54      # Convert to intermediate cartesian representation (N x 3) 
 55      cart_orig = np.empty(shape=(len(sph_coords), 3), dtype=float) 
 56      cart_orig[:, 0] = c(sph_coords[:, 1]) # * s(sph_coords[:, 0]) 
 57      cart_orig[:, 1] = s(sph_coords[:, 1]) # * s(sph_coords[:, 0]) 
 58      cart_orig[:, 0:2] *= s(sph_coords[:, 0])[:, None] 
 59      cart_orig[:, 2] = c(sph_coords[:, 0]) 
 60   
 61      # Rotate by x_i = R_{ij} * x_j 
 62      # NB: dot contracts last dim of A with second-to-last dim of B 
 63      cart_new = np.dot(cart_orig, R) 
 64   
 65      # Extract new spherical coordinates 
 66      sph_new = np.empty(shape=(len(sph_coords), 2), dtype=float) 
 67      sph_new[:, 0] = np.arccos(cart_new[:, 2]) 
 68      sph_new[:, 1] = np.arctan2(cart_new[:, 1], cart_new[:, 0]) 
 69   
 70      return sph_new 
 71   
 72 -def new_z_to_euler(new_z): 
 73      """ 
 74      From the new Z axis expressed in (polar, azimuthal) angles of the 
 75      initial coordinate system, return the (alpha, beta) Euler angles 
 76      that rotate the old Z axis (0, 0) to the new Z axis. 
 77      """ 
 78      return (lal.PI_2 + new_z[1]) % (2 * lal.PI), new_z[0] 
 79   
 80 -def rotate_about_axis(x, axis, ang): 
 81      """ 
 82      Return the result of rotating x about axis by ang (in radians). 
 83      axis must be a unit vector. 
 84      """ 
 85      if abs(np.dot(axis, axis) - 1) >= 1e-6: 
 86          raise ValueError, "axis must be a unit vector" 
 87      if len(x) != 3: 
 88          raise ValueError, "x must be three-dimensional" 
 89      if len(axis) != 3: 
 90          raise ValueError, "axis must be three-dimensional" 
 91      cosa = np.cos(ang) 
 92      sina = np.sin(ang) 
 93   
 94      # Rodrigues' rotation formula 
 95      R = np.array( 
 96          [[cosa+axis[0]*axis[0]*(1.0-cosa),  
 97            axis[0]*axis[1]*(1.0-cosa)-axis[2]*sina, 
 98            axis[0]*axis[2]*(1.0-cosa)+axis[1]*sina], 
 99           [axis[1]*axis[0]*(1.0-cosa)+axis[2]*sina, 
100            cosa+axis[1]*axis[1]*(1.0-cosa), 
101            axis[1]*axis[2]*(1.0-cosa)-axis[0]*sina], 
102           [axis[2]*axis[0]*(1.0-cosa)-axis[1]*sina, 
103            axis[2]*axis[1]*(1.0-cosa)+axis[0]*sina, 
104            cosa+axis[2]*axis[2]*(1.0-cosa)]]) 
105      return np.dot(R, x) 
106  # 
107  # Utilities to find the angle between two points 
108  # 
109   
110 -def _abs_diff(c): 
111      """ 
112      For some angular difference c = |a - b| in radians, find the 
113      magnitude of the difference, taking into account the wrap-around at 2*pi. 
114      """ 
115      c = abs(c) % (2 * lal.PI) 
116      # XXX: numpy 1.3.0 introduces fmin, which is more elegant 
117      return np.where(c < lal.PI, c, 2 * lal.PI - c) 
118   
119 -def _haversine(angle): 
120      return np.sin(angle / 2)**2 
121   
122 -def _archaversine(h): 
123      """ 
124      Compute the inverse of the _haversine function, using clip as protection 
125      for the antipodes. 
126      """ 
127      h = np.clip(h, 0., 1.) 
128      return 2 * np.arcsin(np.sqrt(h)) 
129   
130 -def angle_between_points(a, b): 
131      """ 
132      Find the angle in radians between a and b, each expressed in 
133      (polar, azimuthal) angles in radians. If a and b are Nx2 arrays 
134      of vectors, then return the N pairwise angles between them. 
135   
136      This formula is the law of _haversines, which is derivable from the 
137      spherical law of cosines, but is more numerically stable about a == b. 
138      This technique is slightly unstable for antipodal a and b. 
139      """ 
140      s = np.sin 
141      dtheta, dphi = (a - b).T 
142      h = _haversine(dtheta) + s(a[..., 0]) * s(b[..., 0]) * _haversine(dphi) 
143      result = _abs_diff(_archaversine(h)) 
144      if result.shape == (1,): 
145          return result[0] 
146      else: 
147          return result 
148   
149  # 
150  # Implement the Fisher distribution 
151  # 
152   
153 -def fisher_rvs(mu, kappa, size=1): 
154      """ 
155      Return a random (polar, azimuthal) angle drawn from the Fisher 
156      distribution. Assume that the concentration parameter (kappa) is large 
157      so that we can use a Rayleigh distribution about the north pole and 
158      rotate it to be centered at the (polar, azimuthal) coordinate mu. 
159   
160      Assume kappa = 1 / sigma**2 
161   
162      References: 
163        * http://en.wikipedia.org/wiki/Von_Mises-Fisher_distribution 
164        * http://arxiv.org/pdf/0902.0737v1 (states the Rayleigh limit) 
165      """ 
166      rayleigh_rv = \ 
167          np.array((np.random.rayleigh(scale=1. / np.sqrt(kappa), size=size), 
168                    np.random.uniform(low=0, high=2*lal.PI, size=size)))\ 
169                  .reshape((2, size)).T  # guarantee 2D and transpose 
170      a, b = new_z_to_euler(mu) 
171      return rotate_euler(rayleigh_rv, a, b, 0) 
172   
173 -def fisher_pdf(theta, kappa): 
174      """ 
175      Return the PDF of theta, the opening angle of X with respect to mu where 
176      X is Fisher-distributed about mu. See fisher_rvs for the definition of mu. 
177      """ 
178      return kappa / (2 * np.sinh(kappa)) * np.exp(kappa * np.cos(theta))\ 
179          * np.sin(theta) 
180   
181 -def fisher_cdf(theta, kappa): 
182      return 0.5 * (np.exp(kappa) - np.exp(kappa * np.cos(theta))) / np.sinh(kappa) 
183