pylal.cSqrTutils

1 """ 2 * Copyright (C) 2004, 2005 Cristina V. Torres 3 * 4 * This program is free software; you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation; either version 2 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with with program; see the file COPYING. If not, write to the 16 * Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, 17 * MA 02111-1307 USA 18 """ 19 __author__ = 'Cristina Torres <cristina.torres@ligo.org>' 20 __date__ = '$Date$' 21 __version__ = '' 22 23 24 import os 25 from numpy import ndarray,dtype,NaN,Inf,zeros,mean,sin,\ 26 pi,all,any,array,arange,inner,cross,isnan,arcsin,\ 27 mod 28 from pylal.metaarray import TimeSeries,TimeSeriesList 29 import copy 30 import sys 31 from StringIO import StringIO 32 from commands import getstatusoutput 33 from scipy import interp 34 from followup_utils import connectToFrameData 35 from pylal import metaarray 36 import warnings 37 38 try: 39 import sqlite3 as sqlite 40 except ImportError: 41 import sqlite 42 if os.getenv("DISPLAY") == None: 43 #Non-interactive 44 try: 45 import matplotlib 46 matplotlib.use("Agg") 47 import pylab 48 except Exception, errorInfo: #RuntimeError,ImportError: 49 disableGraphics=True 50 sys.stderr.write("Error trying to import NON-INTERACTIVE pylab!\n") 51 sys.stderr.write("Exception Instance :%s\n"%(str(type(errorInfo)))) 52 sys.stderr.write("Exception Args :%s\n"%(str(errorInfo.args))) 53 sys.stderr.write("Pylab functionality unavailable!\n") 54 else: 55 #Interactive 56 try: 57 import pylab 58 except Exception, errorInfo: #RuntimeError,ImportError: 59 disableGraphics=True 60 sys.stderr.write("Error trying to import INTERACTIVE pylab!\n") 61 sys.stderr.write("Exception Instance :%s\n"%(str(type(errorInfo)))) 62 sys.stderr.write("Exception Args :%s\n"%(str(errorInfo.args))) 63 sys.stderr.write("Pylab functionality unavailable!\n") 64 65 """ 66 This file will hold all the methods required to run the c^2t(Cristina 67 Cesar Triangulation) DetChar tool. 68 """ 69 70 71 # 72 # Iterate in matrix computing row dot vectors 73 #

74 -def CSTdot(A,B):

75 """ 76 Input two matrices Nx3 in size. 77 Calculates the inner product row by row. 78 """ 79 if( A.shape != B.shape) and A.ndim > 2: 80 os.abort() 81 if A.ndim == 2: 82 return array([inner(A[i],B[i]) for i in arange(A.shape[0])]) 83 else: 84 return array(inner(A,B))

85 # 86 # Iterate in matrix computing row cross vectors 87 #

88 -def CSTcross(A,B):

89 """ 90 Input two matrices Nx3 in size. 91 Calculates the cross product row by row. 92 """ 93 if( A.shape != B.shape) and A.ndim > 2: 94 os.abort() 95 if A.ndim == 2: 96 return array([cross(A[i],B[i]) for i in arange(A.shape[0])]) 97 else: 98 return array(cross(A,B))

99 100 # 101 # Intersect two planes 102 #

103 -def intersectPlanes(normA,normB,coordA,coordB,wiggleA=0.0):

104 """ 105 Expects a 3xM matrices of 3D vectors for the plane norms 106 also in addition to that we take vectors from the "origin" 107 (0,0,0) (single 3x1 matrix) to be the coordinate input. The 108 output is a tuple of 3xM arrays parameterized by their index, 109 intersection vector and point on that line. Output is 110 (linePoint,LineVector) ((3xM),(3xM)). The input wiggleA is the 111 wiggle angle in radian to consider two vectors parrallel. 112 """ 113 #If inputs a 'lone' vectors make them a 3x1 matrix! 114 #Normalize the input vectors, but NOT coords! 115 for j in arange(len(normA)): 116 normA[j]=normA[j]/CSTdot(normA[j],normA[j]) 117 normB[j]=normB[j]/CSTdot(normB[j],normB[j]) 118 #Cycle through all the rows of input matrix 119 lineVec=zeros(normA.shape) 120 linePoi=zeros(normA.shape) 121 for i in arange(len(normA)): 122 #Determine if planes parallel? 123 if arcsin(CSTdot(normA[i],normB[i])) >= (1-wiggleA): 124 #Set values to NaN for given 'i' 125 lineVec[i]=NaN 126 linePoi[i]=NaN 127 else: 128 #Find normal vector along line of intersection! 129 tmpVec=CSTcross(normA[i],normB[i]) 130 lineVec[i]=tmpVec/CSTdot(tmpVec,tmpVec) 131 print "Intersecting line vector(unit) %s"%(lineVec) 132 #Find point along line of intersection! 133 if isnan(lineVec[i]).any() or all(lineVec[i]==0.0): 134 print "Vector of plane intersection line is ZERO Vector!" 135 lineVec[i]=NaN 136 linePoi[i]=NaN 137 else: 138 #Determine which zero axis vector points through! 139 zeroIndex=lineVec.argmax() 140 #Map out x,y,z to index of i,j,k with lowest index as 'i' 141 (iX,iY,iZ)=mod(range(zeroIndex,zeroIndex+3),3).tolist() 142 print "Vector elements:%s,%s,%s"%(iX,iY,iZ) 143 #(x,y,z)==(l0,l1,l2) 144 linePoi[i][iZ]=0 145 linePoi[i][iX]=(CSTdot(normB[i],coordB)-(CSTdot(normA[i],coordA)*normB[i][iY])/normA[i][iY])/(normB[i][iX]-(normB[i][iY]*normA[i][iX])/normA[i][iY]) 146 linePoi[i][iY]=(CSTdot(normA[i],coordA)-normA[i][iX]*linePoi[i][iX])/normA[i][iY] 147 return (lineVec,linePoi)

148 149 # 150 # Intersect two lines 151 #

152 -def intersectLines(vecA,vecB,coordA,coordB,wiggleA=0.0):

153 """ 154 Expects inputs of 3xM matrices representing vectors along a line 155 and a vector 3xM of matching points on lines for those vectors. 156 Output is in a 3xM collection of intersecting points. 157 Examples: 158 Input 159 ([[vX,vY,vZ]],[[lX,lY,lZ]]) 160 161 Output 162 [[i,j,k]] 163 A set of lines is considered parallel if angle(rads) between 164 then is less than wiggleA, in this case the output coordinate 165 is (NaN,NaN,NaN). 166 """ 167 #Error check input shapes etc! 168 #Normalize the input vectors, but NOT coords! 169 for j in arange(len(vecA)): 170 vecA[j]=vecA[j]/CSTdot(vecA[j],vecA[j]) 171 vecB[j]=vecB[j]/CSTdot(vecB[j],vecB[j]) 172 #Generate a matrix of triangulated coordinates 173 triPoi=zeros(vecA.shape) 174 for i in arange(len(vecA)): 175 #Check if line vectors are parallel?(vecA&vecB unit vectors) 176 if arcsin(CSTdot(vecA[i],vecB[i]))<=wiggleA: 177 #Set values to NaN for given 'i' 178 triPoi[i]=NaN 179 else: 180 (iX,iY,iZ)=(0,1,2) 181 #Find intersection of the two lines. 182 #P(s)=P_0+s*\vec{P} 183 gammaConst=(\ 184 coordA[i][iY]*vecA[i][iZ] - vecA[i][iZ]*coordB[i][iY] +\ 185 vecA[i][iY]*(coordB[i][iZ]-coordA[i][iZ])\ 186 )/\ 187 (vecA[i][iY]*vecB[i][iZ]+vecA[i][iZ]*vecB[i][iY]) 188 triPoi=coordB+(gammaConst*vecB) 189 return triPoi

190

191 -class cSqrTSensor(object):

192 """ 193 This class is just method globber that works expects as input 194 TimeSeries objects defined by pylal.frutils. This class only 195 works in a cartesian basis! The distance should be measure in units 196 of meters. 197 """

198 - def __init__(self,\ 199 name="Unknown",\ 200 xAxis=None,\ 201 yAxis=None,\ 202 zAxis=None,\ 203 coordinate=(None,None,None),\ 204 swapSign=(False,False,False),\ 205 tType="normal"\ 206 ):

207 """ 208 This object is initialized with 3 data streams (1xM) which should be the 209 intesities of each spatial coordinate seen by the sensor in 210 question. In addition to this information, the cartesian 211 coordinates relative to some origin should be specified. Inside 212 the class the data will be represented as collections of row 213 vectors in a 3xM matrix. The units 214 are not checked between different cSqrTSensor objects. The user 215 should take care to make sure this is sensible! 216 """ 217 self.name=name.strip().upper() 218 self.__tType__=tType.strip().lower() 219 if self.__tType__ != "normal": 220 os.abort() 221 # 222 # Set variable type as REAL8 223 self.dtype='float64' 224 # 225 # Set angular wiggle room to define parralel 226 wiggleAngle=3.0 #Degrees 227 wiggleDistance=25.0 #Meters 228 self.__wiggleA__=None 229 self.__wiggleD__=None 230 self.setParallelMismatchAngle(wiggleAngle) 231 self.setHypersphereMismatch(wiggleDistance) 232 # 233 # Split off the metadata if input is a TimeSeries 234 # 235 if isinstance(xAxis,TimeSeries): 236 self.metaX=xAxis.metadata 237 else: 238 self.metaX=None 239 if isinstance(yAxis,TimeSeries): 240 self.metaY=yAxis.metadata 241 else: 242 self.metaY=None 243 if isinstance(zAxis,TimeSeries): 244 self.metaZ=zAxis.metadata 245 else: 246 self.metaZ=None 247 # 248 # Verify Input data properties if TimeSeries 249 # 250 if all(self.metaX,self.metaY,self.metaZ): 251 if not ( 252 xMetaDict["segments"].intersects_all(yMetaDict["segments"]) \ 253 and \ 254 yMetaDict["segments"].intersects_all(zMetaDict["segments"]) \ 255 and \ 256 zMetaDict["segments"].intersects_all(xMetaDict["segments"]) \ 257 ): 258 raise Exception, \ 259 "Inputs for X,Y,Z data don't cover identical segments!" 260 # 261 # Check sampling rates! Add ability to resample on the fly! 262 # 263 if not(\ 264 xMetaDict["dt"] == \ 265 yMetaDict["dt"] == \ 266 zMetaDict["dt"] \ 267 ): 268 raise Exception, \ 269 "Input data for sampling rates on data is not consistent!" 270 else: 271 sys.stdout.write("Inputs are not TimeSeries ojects can not verify data traits.\n") 272 273 # 274 # Save the coordinate information 275 # (Coordinate 3x1 Array [X,Y,Z] NOT [[X,Y,Z]]) 276 if isinstance(coordinate,ndarray): 277 self.coordinate=coordinate 278 else: 279 self.coordinate=array(coordinate,dtype=self.dtype) 280 # 281 # Save the input data 282 # 283 self.data=zeros((len(xAxis),3),self.dtype) 284 self.data[:,0]=xAxis 285 self.data[:,1]=yAxis 286 self.data[:,2]=zAxis

287 # 288 # End Init 289 # 290

291 - def setParallelMismatchAngle(self,\ 292 wiggleAngle=3.0\ 293 ):

294 """ 295 Enter an angle in degrees to define a tolerance on the 296 parallelness of two lines for determining potential 297 intersection of lines and planes. 298 """ 299 self.__wiggleA__=sin(wiggleAngle/(2.0*pi))#Radians

300

301 - def setHypersphereMismatch(self,\ 302 wiggleDistance=25.0\ 303 ):

304 """ 305 Enter an angle in degrees to define a tolerance on the 306 parallelness of two lines for determining potential 307 intersection of lines and planes. 308 """ 309 self.__wiggleD__=wiggleDistance

310

311 - def triangulate(self,\ 312 otherSensor=[None],\ 313 ):

314 """ 315 Wrapper to pick the correct triangulation machinery. 316 """ 317 if not isinstance(otherSensor,list): 318 otherSensor=[otherSensor] 319 if not all([isinstance(x,cSqrTSensor) for x in otherSensor]): 320 raise Exception, "Expected list of cSqrTSensor objects, got list of other objects!" 321 if self.__tType__ == "normal": 322 return self.__normalTriangulate__(otherSensor) 323 else: 324 os.abort()

325

326 - def __normalTriangulate__(self,\ 327 otherSensor=[None],\ 328 ):

329 """ 330 Assuming the input is normals to 3D planes in 3 space we will 331 intercept the planes resulting in either points in 3D space or 2D 332 lines whichever is the most physically interesting. 333 """ 334 # 335 # Check for name space collisions 336 if self.name in [x.name for x in otherSensor]: 337 raise Exception, "Namespace collision between sensors to triangulate!" 338 # 339 # Perform triangulation 340 # (self,A,B) 341 planeIntersections=dict() 342 pKeyMask="%s_%s" 343 if len(otherSensor) >= 1: 344 #Intersect planes pairwise 345 sensorList=list() 346 sensorList.append(self) 347 sensorList.extend(otherSensor) 348 #Identify uniq pairs 349 for aa,A in enumerate(sensorList): 350 for bb,B in enumerate(sensorList): 351 myKey=pKeyMask%(A.name,B.name) 352 symKey=pKeyMask%(B.name,A.name) 353 if (aa!=bb) and \ 354 not (myKey in planeIntersections.keys()) and \ 355 not (symKey in planeIntersections.keys()): 356 planeIntersections[pKeyMask%(A.name,B.name)]=intersectPlanes(A.data,\ 357 B.data,\ 358 A.coordinate,\ 359 B.coordinate, 360 wiggleA=self.__wiggleA__) 361 else: 362 raise Exception, "Need at least two normals to intersect planes." 363 #In case of only on pair of plane given as inputs 364 if len(planeIntersections.keys()) == 1: 365 myKey=planeIntersections.keys() 366 return planeIntersections[myKey[0]] 367 #If we have a pairs of lines to intersect 368 lineIntersections=dict() 369 for (nameA,lineA) in planeIntersections.iteritems(): 370 for (nameB,lineB) in planeIntersections.iteritems(): 371 lineIntersections[pKeyMask%(nameA,nameB)]=intersectLines(lineA[0],\ 372 lineB[0],\ 373 lineA[1],\ 374 lineB[1],\ 375 wiggleA=self.__wiggleA__) 376 # 377 # If points at time index all inside this hypersphere return mean point! 378 # 379 resultPoint=zeros(lineIntersections[0],dtype=self.dtype) 380 resultVectors=zeros(lineIntersections[0],dtype=self.dtype) 381 #If points close enough return centroid else None 382 for i in arange(len(lineIntersections[0])): 383 tmpCentroid=zeros(lineIntersections[0],dtype=self.dtype) 384 for j in arange(len(lineIntersections)): 385 tmpCentroid+=lineIntersections[j] 386 tmpCentroid=tmpCentroid/float(len(lineIntersections)) 387 resultPoint[i]=tmpCentroid 388 if max([sqrt(CSTdot(x-tmpCentroid,x-tmpCentroid)) \ 389 for x in lineIntersections.iteritems]) > \ 390 self.__wiggleD__: 391 resultPoint[i]=NaN 392 else: 393 resultPoint[i]=tmpCentroid 394 # 395 # Return triangulation coordinates 396 # 397 return (resultVectors,resultPoint)

398

Source Code for Module pylal.cSqrTutils