← Index
Performance Profile   « block view • line view • sub view »
For /wise/base/deliv/dev/bin/getfix
  Run on Thu May 20 15:30:03 2010
Reported on Thu May 20 16:25:29 2010

File/wise/base/deliv/dev/lib/perl/WISE/CoUtils.pm
Statements Executed31121298
Total Time80.8105209982546 seconds

Subroutines — ordered by exclusive time
Calls P F Exclusive
Time		 Inclusive
Time		 Subroutine
59809671161.0436479.71310WISE::CoUtils::angdist
55478110.768070.76807WISE::CoUtils::hdms_dd
27739310.603474.34133WISE::CoUtils::cconv
27739110.271101.03917WISE::CoUtils::to_dd
13869110.140520.24279WISE::CoUtils::rd2radvec
13869110.118912.49860WISE::CoUtils::_equgal
13869110.096142.28587WISE::CoUtils::_fk5fk4
13869110.092310.20001WISE::CoUtils::_equecl
1113.7e-53.7e-5WISE::CoUtils::_fpre
1111.1e-57.7e-5WISE::CoUtils::_eclequ
00000WISE::CoUtils::BEGIN
00000WISE::CoUtils::_eclgal
00000WISE::CoUtils::_fk4fk5
00000WISE::CoUtils::_galecl
00000WISE::CoUtils::_galequ
00000WISE::CoUtils::_trivial
00000WISE::CoUtils::_undefize
00000WISE::CoUtils::autocorrhist
00000WISE::CoUtils::cconvdd
00000WISE::CoUtils::dd_hdms
00000WISE::CoUtils::fit_3param_simple
00000WISE::CoUtils::gcposang
00000WISE::CoUtils::gnomfwd_pure
00000WISE::CoUtils::gnomrev_pure
00000WISE::CoUtils::lonlat2rpa_small
00000WISE::CoUtils::platefit_3param_simple
00000WISE::CoUtils::to_sex
111-152.71548460.76892WISE::CoUtils::posmatch

LineStmts.Exclusive
Time		Avg.Code
1#! /usr/bin/env perl
233.3e-51.1e-5use strict;
# spent 13µs making 1 call to strict::import
335.1e-51.7e-5use warnings;
# spent 36µs making 1 call to warnings::import
4
5package WISE::CoUtils;
6
732.9e-59.7e-6use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $VERSION $R2D $PI $S2);
# spent 98µs making 1 call to vars::import
8
9# Extension modules
10
1137.7e-52.6e-5use Exporter;
# spent 28µs making 1 call to Exporter::import
1211.0e-61.0e-6$VERSION = 1.00;
1311.0e-51.0e-5@ISA = qw(Exporter);
14
1511.0e-61.0e-6@EXPORT = ();
16
1717.0e-67.0e-6@EXPORT_OK = qw(angdist posmatch to_dd to_sex rd2rad
18 rd2radvec lonlat2rpa_small gcposang
19 cconv convdd autocorrhist
20 gnomrev_pure gnomfwd_pure
21 $R2D $PI $S2 $sig2fwhm);
22
2330.000590.00020use WISE::CoUtilsXS;
# spent 154µs making 1 call to Exporter::import
24
2533.7e-51.2e-5use POSIX ('tan','acos', 'asin');
# spent 82µs making 1 call to POSIX::import
26
2730.000540.00018use vars qw /$PI $R2D $S2 $sig2fwhm/;
# spent 58µs making 1 call to vars::import
28
29100$PI = 3.1415926535897931;
3012.0e-62.0e-6$R2D = 180.0 / $PI;
3111.0e-61.0e-6$S2 = sqrt(2.0);
32100$sig2fwhm = 2.354820;
33
3413.0e-63.0e-6my $iam = "$0/CoUtils";
3511.0e-61.0e-6my $err = "*** $iam";
3611.0e-61.0e-6my $warn= "=== $iam";
37
38# CVS revision ID
3911.0e-61.0e-6my $version = '$Id: CoUtils.pm 6953 2010-01-14 01:35:33Z tim $ ';
40
41# Supplement and/or alias WISE::Pop and Astro::Coord conversions.
42# Convert to use in degrees instead of "turns."
43sub _fk4fk5 { return map {$_*360} fk4fk5($_[0]/360,$_[1]/360); }
44138690.099407.2e-6
# spent 2.29s (96.1ms+2.19) within WISE::CoUtils::_fk5fk4 which was called 13869 times, avg 165µs/call: # 13869 times (96.1ms+2.19s) by WISE::CoUtils::_equgal at line 55, avg 165µs/call
sub _fk5fk4 { return map {$_*360} fk5fk4($_[0]/360,$_[1]/360); }
# spent 2.19s making 13869 calls to Astro::Coord::fk5fk4, avg 158µs/call
45# - We assume the user wants J2000 equatorial coords, except when
46# converting explicitly to B1950.
47# - We assume the user wants the J2000 ecliptic.
48# - The galactic routines take/use b1950 ra,dec, so
49# we use b1950 equatorial coord.s as an intermediary.
50sub _galequ { my ($bra,$bdec) = galeqd($_[0],$_[1]);
51 return _fk4fk5($bra,$bdec); }
52sub _galecl { my ($bra,$bdec) = galeqd($_[0],$_[1]);
53 ($bra,$bdec) = _fk4fk5($bra,$bdec);
54 return equecx($bra,$bdec,2,2000); }
55277380.198057.1e-6
# spent 2.50s (119ms+2.38) within WISE::CoUtils::_equgal which was called 13869 times, avg 180µs/call: # 13869 times (119ms+2.38s) by WISE::CoUtils::cconv at line 130, avg 180µs/call
sub _equgal { my ($bra,$bdec) = _fk5fk4($_[0],$_[1]);
# spent 2.29s making 13869 calls to WISE::CoUtils::_fk5fk4, avg 165µs/call
56 return equgad($bra,$bdec); }
# spent 93.8ms making 13869 calls to WISE::CoUtilsXS::equgad, avg 7µs/call
57138690.149631.1e-5
# spent 200ms (92.3+108) within WISE::CoUtils::_equecl which was called 13869 times, avg 14µs/call: # 13869 times (92.3ms+108ms) by WISE::CoUtils::cconv at line 130, avg 14µs/call
sub _equecl { return equecx($_[0],$_[1],2,2000); }
# spent 108ms making 13869 calls to WISE::CoUtilsXS::equecx, avg 8µs/call
5817.0e-57.0e-5
# spent 77µs (11+66) within WISE::CoUtils::_eclequ which was called # once (11µs+66µs) by WISE::CoUtils::cconv at line 130
sub _eclequ { return ecleqx($_[0],$_[1],2,2000); }
# spent 66µs making 1 call to WISE::CoUtilsXS::ecleqx
59sub _eclgal { my ($bra,$bdec) = ecleqx($_[0],$_[1],2,2000);
60 ($bra,$bdec) = _fk5fk4($bra,$bdec);
61 return equgad($bra,$bdec); }
62# For non-converting conversions
63sub _trivial { return @_; }
64
65# Acceptable coordinate system aliases
6612.0e-52.0e-5our %coords = (qw/equ equ
67 eq equ
68 1 equ
69 q equ
70 equatorial equ
71 j2000 equ
72 fk5 equ
73 2000 equ
74 gal gal
75 ga gal
76 g gal
77 2 gal
78 galactic gal
79 ecl ecl
80 ec ecl
81 3 ecl
82 ecliptic ecl
83 b1950 fk4
84 fk4 fk4
85 1950 fk4 /
86 );
87
88# Pair conversion functions
8911.7e-51.7e-5my %conv = (
90 equequ => \&_trivial,
91 equgal => \&_equgal,
92 equecl => \&_equecl,
93 galequ => \&_galequ,
94 galgal => \&_trivial,
95 galecl => \&_galecl,
96 eclequ => \&_eclequ,
97 eclgal => \&_eclgal,
98 eclecl => \&_trivial,
99 fk4equ => \&_fk4fk5,
100 equfk4 => \&_fk5fk4,
101 fk4fk4 => \&_trivial,
102 );
103
10430.001250.00042use Astro::Coord;
# spent 133µs making 1 call to Exporter::import
10530.005860.00195use WISE::CoUtilsXS;
# spent 144µs making 1 call to Exporter::import
106
107
# spent 4.34s (603ms+3.74) within WISE::CoUtils::cconv which was called 27739 times, avg 157µs/call: # 13869 times (271ms+2.98s) at line 838 of /wise/base/deliv/dev/bin/getfix, avg 234µs/call # 13869 times (333ms+758ms) at line 833 of /wise/base/deliv/dev/bin/getfix, avg 79µs/call # once (53µs+237µs) at line 485 of /wise/base/deliv/dev/bin/getfix
sub cconv {
1083883460.690721.8e-6 warn '$0/Cconv: Bad args. cconv($from,$lon,$lat[,$to])'
109 if @_ != 3 && @_ != 4;
110 my $from = shift;
111 my $lon0 = shift;
112 my $lat0 = shift;
113 my $to = shift || 'equ';
114
115 warn("$0/Cconv: Some required args undefined.\n"), return
116 if grep {! defined $_} ($from,$lon0,$lat0);
117 $from = $coords{lc $from}
118 or warn("$0/Cconv: Coord alias '$from' not known.\n"), return;
119 $to = $coords{lc $to}
120 or warn("$0/Cconv: Coord alias '$to' not known.\n"), return;
121 my $hord = $from eq 'equ' ? 'h' : 'd';
122 # Convert from sexigesimal (detected internally in to_dd)
123 my ($lon0d,$lat0d) = to_dd($lon0,$lat0,$hord);
# spent 1.04s making 27739 calls to WISE::CoUtils::to_dd, avg 37µs/call
124 warn("$0/Cconv: Bad position ($lon0,$lat0) passed.\n"), return
125 if ! defined $lon0d || ! defined $lat0d || $lat0d < -90 || $lat0d > 90;
126 warn "$0/Cconv: Unnormalized longitude ($lon0) passed.\n"
127 if $lon0d < 0 || $lon0d > 360;
128 my $conv = $conv{$from.$to}
129 or warn("$0/Cconv: Conversion function '$from$to' not available.\n"), return;
130 return ($conv->($lon0d,$lat0d),$from,$to);
# spent 2.50s making 13869 calls to WISE::CoUtils::_equgal, avg 180µs/call # spent 200ms making 13869 calls to WISE::CoUtils::_equecl, avg 14µs/call # spent 77µs making 1 call to WISE::CoUtils::_eclequ
131}
132
133# Version of above where positions are guaranteed to be decimal degrees
134# and in range
135sub cconvdd {
136 my $from = shift;
137 my $lon0 = shift;
138 my $lat0 = shift;
139 my $to = shift || 'equ';
140
141 $from = $coords{lc $from}
142 or warn("$0/Cconv: Coord alias '$from' not known"), return;
143 $to = $coords{lc $to}
144 or warn("$0/Cconv: Coord alias '$to' not known"), return;
145 my $conv = $conv{$from.$to}
146 or warn("$0/Cconv: Conversion function '$from$to' not available"), return;
147 return ($conv->($lon0,$lat0),$from,$to);
148}
149
150# Pure perl gnomonic projections
151sub gnomrev_pure {
152 my ($xi,$eta,$ra0,$dec0,$crota,$cdelt,$cenpix) = @_;
153 my $x = ($xi - $cenpix)*$cdelt/$R2D;
154 my $y = ($eta - $cenpix)*$cdelt/$R2D;
155 my $crotar = $crota/$R2D;
156 my $xp = $x*cos($crotar) + $y*sin($crotar);
157 my $yp = -$x*sin($crotar) + $y*cos($crotar);
158 # delta = atan(sqrt($xp*$xp + $yp*$yp))
159 my $rsq = $xp*$xp + $yp*$yp;
160 my $r = sqrt($rsq);
161 my $cosdelta = $rsq > 1e-32 ? 1. / sqrt($rsq + 1) : 1.0;
162 my $sindelta = sqrt(1. - $cosdelta*$cosdelta);
163 # beta = atan2(-$x,$y)
164 my $cosbeta = $r>0 ? -$yp/$r : 1;
165 my $sinbeta = $r>0 ? -$xp/$r : 0;
166 my $cosdec0 = cos($dec0/$R2D);
167 my $sindec0 = sin($dec0/$R2D);
168 my $xx = $sindec0*$sindelta*$cosbeta + $cosdec0*$cosdelta;
169 my $yy = $sindelta*$sinbeta;
170 my $dec = asin($sindec0*$cosdelta - $cosdec0*$sindelta*$cosbeta)*$R2D;
171 my $ra = $ra0 + atan2($yy,$xx)*$R2D;
172 if($ra >= 360) { $ra -= 360 ; }
173 if($ra < 0) { $ra += 360 ; }
174 #printf("--- rev: returning %f,%f.\n",$ra,$dec);
175 return ($ra,$dec);
176}
177
178sub gnomfwd_pure {
179 my ($ra, $dec, $ra0, $dec0, $crota, $cdelt, $cenpix) = @_;
180 my $cosdec0 = cos($dec0/$R2D);
181 my $sindec0 = sin($dec0/$R2D);
182 my $drar = ($ra0-$ra)/$R2D;
183 my $decr = $dec/$R2D;
184 my $crotar = $crota/$R2D;
185 my $a = cos($decr) * cos($drar);
186 my $f = 1.0/($sindec0 * sin($decr) + $a*$cosdec0);
187 my $xp = $f * cos($decr) * sin($drar);
188 my $yp = $f * ($cosdec0 * sin($decr) - $a*$sindec0);
189 my $x = $xp*cos(-$crotar) + $yp*sin(-$crotar);
190 my $y = -$xp*sin(-$crotar) + $yp*cos(-$crotar);
191 my $xi = $x*$R2D/$cdelt + $cenpix;
192 my $eta = $y*$R2D/$cdelt + $cenpix;
193 return ($xi,$eta);
194}
195
196# Angular distance between two lon-lat pairs. Not great for really small
197# differences.
198
# spent 79.7s (61.0+18.7) within WISE::CoUtils::angdist which was called 5980967 times, avg 13µs/call: # 5980967 times (61.0s+18.7s) by WISE::CHealPix::pix_withinr_ring at line 89 of /wise/base/deliv/dev/lib/perl/x86_64-linux-thread-multi/WISE/CHealPix.pm, avg 13µs/call
sub angdist {
1992990483566.024292.2e-6 my ($x1,$y1,$x2,$y2) = map {$_/$R2D} @_;
200 my $arg = sin($y1)*sin($y2) + cos($y1)*cos($y2)*cos($x1-$x2);
201 $arg = 1 if $arg > 1;
202 $arg = -1 if $arg < -1;
203 return acos($arg)*$R2D;
# spent 18.7s making 5980967 calls to POSIX::acos, avg 3µs/call
204}
205
206
# spent 243ms (141+102) within WISE::CoUtils::rd2radvec which was called 13869 times, avg 18µs/call: # 13869 times (141ms+102ms) by WISE::CoUtils::posmatch at line 767, avg 18µs/call
sub rd2radvec {
207693450.206993.0e-6 my ($r,$dx,$dy) = radifs($_[2],$_[3],$_[0],$_[1]); # From CoUtilsXS
# spent 102ms making 13869 calls to WISE::CoUtilsXS::radifs, avg 7µs/call
208 my $pa = atan2($dy,-$dx)*$R2D-90;
209 $pa += 360 if $pa < 0;
210 $pa -= 360 if $pa > 360;
211 return wantarray ? ($r/3600,$pa,$dx/3600,$dy/3600) : $r/3600;
212}
213
214sub lonlat2rpa_small {
215 my ($ra,$dec,$ra0,$dec0) = @_;
216 my ($x,$y) = gnomfwd_pure($ra,$dec,$ra0,$dec0,0,1/3600,0);
217 my ($r,$pa) = xy2rve(-$x,$y,0,0);
218 return wantarray ? ($r/3600,$pa) : $r/3600;
219}
220
221# Formal position angle using angle between great circle intersecting
222# the two positions and the meridian.
223# Cribbed from IDL posang routine.
224sub gcposang {
225 # The order of the args is switched to match the IDL formula
226 my ($rarad2,$decrad2,$rarad1,$decrad1) = map { $_/$R2D } @_;
227 my $radif = $rarad2-$rarad1;
228 my $angle = atan2(sin($radif),
229 cos($decrad1)*tan($decrad2)-sin($decrad1)*cos($radif))
230 * $R2D;
231 $angle += 360 if $angle < 0;
232 return $angle;
233}
234
235sub fit_3param_simple {
236 my ($ra1, $dec1) = (shift, shift);
237 my ($ra0, $dec0) = (shift, shift); # Truth
238 my ($racen, $deccen) = (shift, shift); # Proj center
239 my $n = @$ra1;
240
241 my (@x1, @y1);
242 for (0..$n-1) {
243 my ($x,$y) = gnomfwd_pure($ra1->[$_], $dec1->[$_],
244 $racen, $deccen, 0, 1/3600., 0);
245 push @x1, $x;
246 push @y1, $y;
247 }
248
249 my (@x0, @y0);
250 for (0..$n-1) {
251 my ($x,$y) = gnomfwd_pure($ra0->[$_], $dec0->[$_],
252 $racen, $deccen, 0, 1/3600., 0);
253 push @x0, $x;
254 push @y0, $y;
255 }
256
257 my $res = platefit_3param_simple(\@x1, \@y1, \@x0, \@y0);
258
259 return $res;
260}
261
262sub platefit_3param_simple {
263 my ($x_in, $y_in) = (shift, shift);
264 my ($xi, $eta) = (shift, shift);
265 my $cdelt = shift;
266
267 require PDL::Lite;
268
269 if(ref($x_in) =~ /array/i) {
270 $x_in = PDL->new($x_in);
271 $y_in = PDL->new($y_in);
272 $xi = PDL->new($xi);
273 $eta = PDL->new($eta);
274 }
275
276 return if $x_in->nelem < 2;
277
278 my ($x_scale,$y_scale) = (1,1);
279 if($cdelt) {
280 $x_scale = abs($cdelt)*3600;
281 $y_scale = abs($cdelt)*3600;
282 }
283
284 # Rescale x_in and y_in and populate x and y
285 my $x = $x_scale*$x_in;
286 my $y = $y_scale*$y_in;
287
288 # Compute the centroids, which will later be used to compute offsets
289 my $x0 = $x->avg;
290 my $y0 = $y->avg;
291
292 my $xi0 = $xi->avg;
293 my $eta0 = $eta->avg;
294
295 # Compute the sums for the rotation angle.
296 my $x_xi = ( ($x-$x0) * ($xi -$xi0) )->sum;
297 my $x_eta = ( ($x-$x0) * ($eta-$eta0) )->sum;
298 my $y_xi = ( ($y-$y0) * ($xi -$xi0) )->sum;
299 my $y_eta = ( ($y-$y0) * ($eta-$eta0) )->sum;
300
301 # Determinant
302 my $det = $x_xi*$y_eta - $x_eta*$y_xi;
303
304 my ($num,$denom);
305 if($det < 0) {
306 $num = $y_xi + $x_eta;
307 $denom = -$x_xi + $y_eta;
308 } else {
309 $num = $y_xi - $x_eta;
310 $denom = $x_xi + $y_eta;
311 }
312
313 # Rotation angle
314 my $theta = atan2( $num, $denom );
315 $theta += 2*$PI if $theta < 0;
316
317 # Package up the fit parameters
318
319 my $ctheta = cos( $theta );
320 my $stheta = sin( $theta );
321
322 my ($ctheta_x,$stheta_y);
323 if($det < 0) {
324 $ctheta_x = -$ctheta;
325 $stheta_y = -$stheta;
326 } else {
327 $ctheta_x = $ctheta;
328 $stheta_y = $stheta;
329 }
330
331 my $stheta_x = $stheta;
332 my $ctheta_y = $ctheta;
333
334 # 3-param solution
335 my $rot = $theta*$R2D;
336 my $xoff = $xi0 - ( $ctheta_x*$x0 + $stheta_x*$y0);
337 my $yoff = $eta0 - (-$stheta_y*$x0 + $ctheta_y*$y0);
338
339 # 6-param solution
340 my @params;
341 $params[0] = $xoff; # a
342 $params[1] = $ctheta_x*$x_scale; # b
343 $params[2] = $stheta_x*$y_scale; # c
344 $params[3] = $yoff; # d
345 $params[4] = -$stheta_y*$x_scale; # e
346 $params[5] = $ctheta_y*$y_scale; # f
347
348 return wantarray ? ($xoff, $yoff, $rot, @params) :
349 {xoff=>$xoff, yoff=>$yoff, rot=>$rot,
350 params=>\@params};
351
352}
353
354sub autocorrhist {
355 my $rows = shift; # Table to autocorrelate
356 my $rmax = shift || 10*60; # Search radius
357 my $opts = shift || {}; # Options my $verbose = $opts->{verbose};
358 my $verbose = $opts->{verbose};
359 my $xy = $opts->{xy};
360 my $debug = $opts->{debug} || "";
361 my $loncol = $opts->{loncol}
362 || $xy ? 'X_IMAGE' : 'X_WORLD'; # Name of the longitude column
363 my $latcol = $opts->{latcol}
364 || $xy ? 'Y_IMAGE' : 'Y_WORLD'; # Name of the latitude column
365 my $magcol = $opts->{magcol}
366 || 'MAG_AUTO'; # Name of the magnitude column
367 my $rbin = $opts->{rbinsz} || 5; # Arc-secs
368 my $pabin = defined $opts->{pabinsz} ? $opts->{pabinsz} : 4; # Degrees
369 my $xybin = $opts->{xybinsz} || 5;
370 my $top = $opts->{top} || 0;
371 my $fits = $opts->{fits}; # FITS input file handle
372 my $norm = $opts->{norm};
373 my $sig = $opts->{sig};
374 my $x0 = $opts->{x0};
375 my $y0 = $opts->{y0};
376 my $fovr = $opts->{fovr};
377 my $pixsz = $opts->{pixsz};
378 my $asymmetric = $opts->{asymmetric}; # Asymmetric matching
379 my $hist2d;
380 my $havecenter;
381 my $top2 = 0;
382
383 if(ref $top) {
384 if(@$top == 1) { $top = $top->[0]; }
385 else {
386 ($top,$top2) = sort {$a<=>$b} @$top
387 if @$top > 1 && $top->[0];
388 }
389 }
390
391 $asymmetric = 0 if ! $top || $top < 0;
392 $asymmetric = 1 if $top2 && $top2 > 0;
393
394 print "Asymmetric matching ($top, $top2) for autocorrelation ".
395 "will be done.\n"
396 if $asymmetric && $verbose;
397
398 if(! ref $rows) {
399 require WISE::FITSIO;
400 $fits = WISE::FITSIO->new($rows,{mode=>'r'})
401 or die;
402 print "Reading from $rows ...\n" if $verbose;
403 $rows = $fits->readcols([$loncol,$latcol,($magcol?($magcol):())],
404 {hdunum=>2})
405 or return;
406 print "... read ".$fits->nread." rows.\n" if $verbose;
407 }
408
409 if($fits && $fovr && $fovr > 0 && ! defined $x0) {
410 $pixsz ||= $fits->key("CDELT2", {hdunum=>1})*3600
411 || $fits->key("PXSCAL2",{hdunum=>1});
412 ($x0,$y0) = $fits->key(["CRVAL1","CRVAL2"],{hdunum=>1});
413 if($xy) {
414 ($x0,$y0) = $fits->cel2pix($x0,$y0);
415 }
416 print "Read FOV center as $x0, $y0 .\n"
417 if $verbose && defined $x0 && defined $y0;
418 }
419
420 ($x0,$y0) = (undef,undef) if defined $x0 && ($x0 < -990 || $y0 < -990);
421
422 my $nrows = @{$rows->{$loncol}};
423
424 if(defined $fovr && $fovr>0 && defined $x0 && defined $y0) {
425 my %newrows;
426 print "Restricting sources to R(FOV) = $fovr\" ...\n" if $verbose;
427 for my $i (0..$nrows-1) {
428 my $d = ($xy
429 ? sqrt(($rows->{$loncol}[$i]-$x0)**2 +
430 ($rows->{$latcol}[$i]-$y0)**2)*$pixsz
431 : angdist($rows->{$loncol}[$i],$rows->{$latcol}[$i],
432 $x0, $y0) * 3600);
433 next if $d > $fovr;
434 push @{$newrows{$_}}, $rows->{$_}[$i] for keys %$rows;
435 }
436 print "... found ".@{$newrows{$loncol}}." rows.\n" if $verbose;
437 $rows = \%newrows;
438 $havecenter = 1;
439 }
440
441 my $mchto = $rows;
442 my $mchfrom = $rows;
443
444 if($magcol && $top && $top <= $#{$rows->{$loncol}}) {
445 my @order;
446 @order = (sort {$rows->{$magcol}[$a]<=>$rows->{$magcol}[$b]}
447 0..$#{$rows->{$loncol}});
448 if($top2 && $top2 <= $#{$rows->{$loncol}}) {
449 print "Restricting match-to group to brightest $top2 sources ...\n"
450 if $verbose;
451 my %newrows;
452 $newrows{$_} = [@{$rows->{$_}}[@order[0..$top2-1]]]
453 for keys %$rows;
454 $mchto = \%newrows;
455 }
456 print "Restricting matching to brightest $top sources ...\n"
457 if $verbose;
458 my %newrows;
459 $newrows{$_} = [@{$rows->{$_}}[@order[0..$top-1]]]
460 for keys %$rows;
461 $mchfrom = \%newrows;
462 $mchto = $mchfrom if ! $asymmetric;
463 }
464
465 # Initialize the histograms. The pre-filled size is important since
466 # only pre-allocated bins are filled.
467
468 my (@rhist, @rnorm, @pahist, @r, @pa, @errnorm, @errpa, @hist2d,
469 @xyhist, @errxy, $nrbins, $npabins, $nxybins);
470
471 $nrbins = int($rmax/$rbin) - 1;
472 @rhist[0.. $nrbins-1] = (0) x $nrbins;
473 @r = map { $_*$rbin } 0..$#rhist;
474
475 if($pabin > 0) {
476 $npabins = int(360/$pabin) - 1;
477 @pahist[0..$npabins-1] = (0) x $npabins;
478 @pa = map { $_*$pabin } 0..$#pahist;
479 } else {
480 $npabins = 0;
481 }
482
483 if($xybin > 0) {
484 $nxybins = 2*int($rmax/$xybin) - 1;
485 @xyhist[0..$nxybins-1] = map { [(0) x $nxybins] } 0..$nxybins-1;
486 }
487
488 # Do the matching and fill the histograms
489
490 my $state = {
491 index => 1,
492 loncol0 => $loncol,
493 latcol0 => $latcol,
494 magcol0 => $magcol,
495 loncol => $loncol,
496 latcol => $latcol,
497 magcol => $magcol,
498 noself => 1,
499 nearest => 0,
500 rhist => \@rhist, rbinsz => $rbin,
501 pahist => \@pahist, pabinsz => $pabin,
502 xyhist => \@xyhist, xybinsz => $xybin,
503 histonly => 1,
504 xy => $xy,
505 pixsz => $pixsz,
506 %$opts
507 };
508
509 my $res = posmatch($mchfrom,$mchto,$xy?$rmax/$pixsz:$rmax,$state);
510 my $n = $res->{nmatched};
511
512 # Fix up the returned histos in various ways
513
514 my $av = $n/($PI*($rmax/3600)**2);
515 for my $i (0..$#rhist) {
516 my $r = ($i+1)*$rbin/3600;
517 # Annular bin area. Corrected for ave. fraction off the fov.
518 my $a = $PI*($r**2 - ($r-$rbin/3600)**2);
519 $rnorm[$i] = $rhist[$i] / $a;
520 $errnorm[$i] = sqrt($rhist[$i])/$a;
521 $rnorm[$i] = $errnorm[$i]>0 ? ($rnorm[$i]-$av)/$errnorm[$i] : 0
522 if $sig;
523 }
524 @rhist = @rnorm if $norm;
525
526 for my $i (0..$#pahist) {
527 $errpa[$i] = sqrt($pahist[$i]);
528 }
529
530 if($xybin > 0) {
531 for my $i (0..$nxybins-1) {
532 for my $j (0..$nxybins-1) {
533 $errxy[$i][$j] = sqrt($xyhist[$i][$j]);
534 }
535 }
536 }
537
538 print "Found $n matches for autocorrelation.\n" if $verbose;
539
540 return wantarray ? @rhist : {
541 n=>$n,
542 nrbins=>$nrbins, npabins=>$npabins,
543 rbinsz=>$rbin, pabinsz=>$pabin,
544 ave=>$av, r=>\@r, pa=>\@pa,
545 rhist=>\@rhist, pahist=>\@pahist,
546 rnorm=>\@rnorm, errnorm=>\@errnorm,
547 errpa=>\@errpa,
548 xybinsz=>$xybin,nxybins=>$nxybins,
549 xyhist=>\@xyhist,errxy=>\@errxy,
550 };
551}
552
553# Do position matching betwee two 'tables' represented in memory as a
554# hash of arrays of columns.
555
# spent 461s (-152715476400ns+613) within WISE::CoUtils::posmatch which was called # once (-152715476400ns+613s) at line 936 of /wise/base/deliv/dev/bin/getfix
sub posmatch {
556751.193620.01591 my $rows = shift; # Table to match from.
557 my $rows0 = shift; # Table to match to.
558 my $rad = shift; # Search radius
559 my $opts = shift || {}; # Options
560 my ($loncol,$latcol,$loncol0,$latcol0,$magcol,$magcol0);
561 my ($hdunum,$index,$nsides,$nearest,$xy,$magtol);
562 my ($nrows,$nrows0,$debug,$verbose,$noself,$pixsz,$histonly);
563 my ($rhist,$rbin,$rmax,$pahist,$pabin,$pamax,$xyhist,$xybin,$xymax);
564
565 if(! $rows) {
566 die "*** $0/POSMCHN: No match-from data rows supplied.\n";
567 }
568
569 if(! $rows0) {
570 die "*** $0/POSMCHN: No match-to data rows supplied.\n";
571 }
572
573 require WISE::CHealPix;
574 WISE::CHealPix->import(qw(ang2pix_ring neighbors_ring pix_withinr_ring
# spent 172µs making 1 call to Exporter::import
575 nsides2side_deg side2nsides_deg));
576
577 $loncol = $opts->{loncol} || 'ra'; # Name of the longitude "column"
578 $latcol = $opts->{latcol} || 'dec'; # Name of the lat "column"
579 $magcol = $opts->{magcol} || 'mag'; # Name of the magnitude column
580 $loncol0 = $opts->{loncol0} || $loncol;
581 $latcol0 = $opts->{latcol0} || $latcol;
582 $magcol0 = $opts->{magcol0} || $magcol;
583 $pixsz = $opts->{pixsz} || 1; # Pixel size for xy mode (for hists)
584 $debug = $opts->{debug} || 0;
585 $hdunum = $opts->{hdunum} || 2;
586 $nsides = $opts->{nsides}; # HealPix nsides param for indexing
587 $verbose = $opts->{verbose};
588 $nearest = $opts->{nearest}; # Return only a single, nearest match
589 $magtol = $opts->{magtol}; # Mag tolerance
590 $noself = $opts->{noself}; # Disallow identical position matches
591 $index = $opts->{index}; # Saved index
592 $rhist = $opts->{rhist};
593 $rmax = $opts->{rmax};
594 $rbin = $opts->{rbinsz};
595 $pahist = $opts->{pahist};
596 $pabin = $opts->{pabinsz};
597 $xyhist = $opts->{xyhist};
598 $xybin = $opts->{xybinsz};
599 $histonly= $opts->{histonly};
600 my $hist = (($rhist && ref($rhist) =~/array/i && $rbin ) ||
601 ($pahist && ref($pahist)=~/array/i && $pabin) ||
602 ($xyhist && ref($xyhist)=~/array/i && $xybin));
603 my $fovr = $xybin*@$xyhist/2 if $hist && $xyhist;
604
605 ++$verbose if $debug;
606
607 die "*** $0/POSMCHN: Radius '"._undefize($rad)."' illegal.\n"
608 if ! defined $rad || ref $rad || $rad<0;
609
610 die "*** $0/POSMCHN: Histogram spec. error; /".
611 join("/", _undefize($rhist,$rbin,$pahist,$pabin,$xyhist,$xybin))."/\n"
612 if ! $hist && ($rhist || $pahist || $xyhist);
613
614 die "*** $0/POSMCHN: lon/latcol/0 name(s) not specified.\n"
615 if grep {! $_} ($loncol,$latcol,$loncol0,$latcol0);
616
617 die "*** $0/POSMCHN: Magtol defined but magcol or magcol0 not specified.\n"
618 if $magtol && (! $magcol || ! $magcol0);
619
620 my %results;
621 my @matched;
622 $results{nmatched} = 0;
623 $results{nmatches} = 0;
624 $results{matched} = \@matched;
625
626 if(! ref $rows) {
627 # Filename was passed, so build rows
628 print "Posmatch: Reading '$rows' ...\n" if $verbose;
629 my $fits = WISE::FITSIO->new($rows,{verbose=>$verbose>1}) or return;
630 $rows = $fits->readcols([$loncol, $latcol, $magtol?($magcol):()],
631 {hdunum=>$hdunum}) or return;
632 }
633
634 if(! ref $rows0) {
635 # Filename was passed, so build rows
636 print "Posmatch: Reading '$rows0' ...\n" if $verbose;
637 my $fits = WISE::FITSIO->new($rows0,{verbose=>$verbose>1}) or return;
638 $rows0 = $fits->readcols([$loncol0, $latcol0, $magtol?($magcol0):()],
639 {hdunum=>$hdunum}) or return;
640 }
641
642 { my @missing =
643 (grep {! $rows->{$_} } ($loncol, $latcol, $magtol?($magcol) :()),
64428.0e-64.0e-6 grep {! $rows0->{$_}} ($loncol0,$latcol0,$magtol?($magcol0):()) );
645 die "*** $0/POSMCHN: Required columns not in data; @missing.\n"
646 if @missing;
647 }
648
649 if(ref($rows) eq 'ARRAY') {
650 # We have a single position as [$ra, $dec]. Turn it into a fake
651 # row structure of the appropriate type.
652 $rows = { $loncol=>[$rows->[0]], $latcol=>[$rows->[1]] };
653 }
654
655 $nrows0 = @{$rows0->{$latcol0}} ;
656 $nrows = @{$rows ->{$latcol} } ;
657
658 if($nrows0 < 1) { return 0; }
659
660 $index = 0 if $nrows0 < 10;
661
66252.0e-54.0e-6 if($index && (! $nsides || $nsides < 0)) {
663 my $side = $rad*2;
664 if($side < 10) { $side = 10; }
665 if($side > 36000) { $side = 36000; }
666 $nsides = side2nsides_deg($side/3600);
# spent 12µs making 1 call to WISE::CHealPix::side2nsides_deg
667 $nsides = 8192 if $nsides>8192;
668 }
669
670 if($hist) {
671 if($rhist) {
672 if(@$rhist) { # Keep pre-allocated size
673 $rmax = $rbin*@$rhist;
674 } else {
675 $rmax ||= $rad;
676 @$rhist = (0) x int($rmax/$rbin);
677 }
678 }
679 if($pahist) {
680 if(@$pahist) {
681 $pamax = $pabin*@$pahist;
682 } else {
683 $pamax ||= 360;
684 @$pahist = (0) x int($pamax/$pabin);
685 }
686 }
687 if($xyhist) {
688 if(@$xyhist) {
689 die "*** $0/POSMCHN: xyhist misconfigured.\n"
690 if ! ref($xyhist->[0])!~/array/i ||
691 grep {@$_ != @$xyhist} @$xyhist;
692 $xymax = $xybin*@$xyhist/2;
693 } else {
694 $xymax ||= $rad;
695 my $nxy = 2*int($xymax/$xybin);
696 $xymax = $xybin*$nxy/2;
697 @$xyhist = map {[(0) x $nxy]} (1..$nxy);
698 }
699 }
700 }
701
702 my $lons0 = $rows0->{$loncol0};
703 my $lats0 = $rows0->{$latcol0};
704 my $lons = $rows ->{$loncol};
705 my $lats = $rows ->{$latcol};
706
707 # If index is defined but is not a reference or is an empty reference,
708 # create a spatial index on the input catalog
70964.1e-56.8e-6 if($index && (! ref($index) || ! keys(%$index) )) {
710 my $side = nsides2side_deg($nsides)*3600;
# spent 23µs making 1 call to WISE::CHealPix::nsides2side_deg
711 print "Posmatch: Indexing with nsides=$nsides, pixel=$side\" ...\n"
712 if $verbose;
713 $index = {} if ! ref($index);
714 %$index = ();
715 # Index rows0
716 for my $ix (0..($nrows0-1)) {
717970832.503192.6e-5 next if ! defined $lons0->[$ix] || ! defined $lats0->[$ix];
718 my $pix = ang2pix_ring($nsides,$lons0->[$ix],$lats0->[$ix]);
# spent 260ms making 13869 calls to WISE::CHealPix::ang2pix_ring, avg 19µs/call
719 my @pixs;
720138692.944510.00021 if($rad <= $side) {
721 @pixs= ($pix,neighbors_ring($nsides,$pix));
722 } else {
723 @pixs= pix_withinr_ring($nsides,$pix,$rad);
# spent 613s making 13869 calls to WISE::CHealPix::pix_withinr_ring, avg 44.2ms/call
724 }
725 print "--- $ix: $lons0->[$ix],$lats0->[$ix] = $pix,@pixs\n"
726 if $debug > 2;
727 push @{ $index->{$_} }, $ix for (@pixs);
728 }
72915.576785.57678 print "Posmatch: ... created ".(keys %$index)." indicies.\n"
730 if $verbose;
731 }
732
733 # If the radius is impossible, just index and return
734 if($rad < 0) { return 0; }
735
736 my @search;
737 @search = (0..$nrows0-1) if ! $index;
738 my $nmatches = 0;
739 my $nmatched = 0;
740
741 print "Posmatch: Comparing $nrows rows to $nrows0 rows ".
742 "with radius $rad\" ...\n"
743 if $verbose;
744
745 SEARCH: for my $ix (0..($nrows-1)) {
746
747100.000555.5e-5 my ($lon,$lat) = ($lons->[$ix], $lats->[$ix]);
748 next if ! defined $lon || ! defined $lat;
749 my @matches;
750 print "\t--- Looking for $lon, $lat ...\n"
751 if $debug > 1;
75230.005870.00196 if($index) {
753 # If there's an index, use it to get a list of rows to check
754 my $pix = ang2pix_ring($nsides,$lon,$lat);
# spent 21µs making 1 call to WISE::CHealPix::ang2pix_ring
755 @search = @{ $index->{$pix} || [] };
756 print "\t\t--- pix=$pix, search=@search\n"
757 if $debug > 1;
758 }
759
760 print "\t\t--- Searching ".@search." rows.\n"
761 if $debug > 1;
762
763 for my $ix0 (@search) {
7641483620.246011.7e-6 my ($lon0,$lat0) = ($lons0->[$ix0], $lats0->[$ix0]);
765 next if ! defined $lon0 || ! defined $lat0;
766 next if $noself && $lon==$lon0 && $lat==$lat0;
767 my($r,$pa) = ( $xy
# spent 243ms making 13869 calls to WISE::CoUtils::rd2radvec, avg 18µs/call
768 ? (sqrt(($lon-$lon0)**2 + ($lat-$lat0)**2),
769 atan2(($lat-$lat0)/$R2D,($lon-$lon0)/$R2D)*$R2D
770 )
771 : rd2radvec($lon,$lat,$lon0,$lat0)
772 ) ;
773 $pa += 360 if $pa < 0;
774 $r *= 3600 if ! $xy;
775 print "\t\t\t--- Comparing to $lon0, $lat0 yields r,pa=$r\",$pa.\n"
776 if $debug > 2;
777 next if $r > $rad;
778 # Mag tolerance check.
779 next if $magtol &&
780 abs($rows0->{$magcol0}[$ix0] - $rows->{$magcol}[$ix])
781 > $magtol;
782 $results{counts2}{$ix0}++;
783 if(! $nearest) {
784 push @matches, {ix=>$ix0,r=>$r,pa=>$pa};
785 } elsif(! @matches || $r < $matches[0]{r}) {
786 $matches[0] = {ix=>$ix0,r=>$r,pa=>$pa};
787 }
788 }
789 $results{counts1}{$ix} = @matches;
790 print "\t\t--- MATCHED ".@matches." rows to #$ix=$lon,$lat".
791 (@matches
792 ?": ".join(";",map {"ix0=$_->{ix},r=$_->{r},".
793 "pos=$lons0->[$_->{ix}],$lats0->[$_->{ix}]"}
794 @matches)
795 :" .")."\n"
796 if $debug > 1 || ($debug && @matches);
797
798 # Record results
79950.001560.00031 if(@matches) {
800
801 $nmatches += @matches;
802 push @matched, $ix;
803 ++$nmatched;
804 $results{matches}{$ix} = [@matches] if ! $histonly;
805
806 if($hist) {
807 my ($bin,$xbin,$ybin);
808
809 for my $mch (@matches) {
810 my ($i,$r,$pa) = @{$mch}{qw/ix r pa/};
811 my $bin;
812 if($rhist && $r < $rmax) {
813 $bin = int($r/$rbin);
814 $rhist->[$bin]++;
815 }
816 if($pahist && $pa < $pamax) {
817 $bin = int($pa/$pabin);
818 $pahist->[$bin]++;
819 }
820 if($xyhist) {
821 my ($x,$y);
822 if(! $xy) {
823 $x = $r*cos($pa/$R2D) + $xymax;
824 $y = -$r*sin($pa/$R2D) + $xymax;
825 } else {
826 $x = $pixsz*$x + $xymax;
827 $y = $pixsz*$y + $xymax;
828 }
829 if($x>=0 && $x<2*$xymax && $y>=0 && $y<2*$xymax) {
830 $xbin = int($x/$xybin);
831 $ybin = int($y/$xybin);
832 $xyhist->[$ybin][$xbin]++;
833 }
834 }
835 } # for
836
837 } # $hist
838
839 } # @matches
840
841 } # SEARCH
842
843 $results{nmatches} = $nmatches;
844 $results{nmatched} = $nmatched;
845 $results{matched} = \@matched;
846
847 if($hist) {
848 if($rhist && @$rhist) { $_||=0 for @$rhist; }
849 if($pahist && @$pahist) { $_||=0 for @$pahist; }
850 if($xyhist && @$xyhist) { $_||=0 for map { (@$_) } @$xyhist; }
851 }
852
853 print "Posmatch: Matched $results{nmatched} source rows to ".
854 "$results{nmatches} target rows.\n"
855 if $verbose;
856
857 return wantarray ? @matched : \%results;
858}
859
860# Convert a sexagesimal (lon,lat) pair representation to decimal degrees
861
# spent 1.04s (271ms+768ms) within WISE::CoUtils::to_dd which was called 27739 times, avg 37µs/call: # 27739 times (271ms+768ms) by WISE::CoUtils::cconv at line 123, avg 37µs/call
sub to_dd {
862 # v => is the lon in hours or degrees?
863554780.318945.7e-6 my ($lon,$lat,$hord) = @_;
864
865 (&hdms_dd($lon,$hord),&hdms_dd($lat,'d'));
# spent 768ms making 55478 calls to WISE::CoUtils::hdms_dd, avg 14µs/call
866}
867
86811.2e-51.2e-5my $FPre = _fpre();
# spent 37µs making 1 call to WISE::CoUtils::_fpre
869
870# Convert a single lon or lat in sexagesimal hours or degrees
871
# spent 768ms within WISE::CoUtils::hdms_dd which was called 55478 times, avg 14µs/call: # 55478 times (768ms+0) by WISE::CoUtils::to_dd at line 865, avg 14µs/call
sub hdms_dd {
8723883460.641671.7e-6 my ($x,$hord) = @_;
873 my ($hd,$m,$s,$pord,$hord2) = ();
874
875 $hord = lc ($hord || "d");
876
877 # Expand compact hhmmss.ss or hh:mm:ss.ss format
878 $x =~ s/^\s*([-+]?)(\d\d?):(\d\d?):(\d\d?(?:\.\d*)?)/$1$2$hord$3m$4s/;
879 $x =~ s/^\s*([-+]?)(\d\d)(\d\d)(\d\d(?:\.\d*)?)/$1$2$hord$3m$4s/;
880
881 if ($x =~ /^$FPre$/) { # F.P. dec. deg.; nothing to do
882 } elsif ($x =~ /[hdms]/) { # Sexigesimal
883 # Confirm correct hms or dms format
884 if($x =~ m'^\s*
885 ([-+])?\s*(?:(\d{1,3}(?:\.\d*)?)([hd]))?\s*
886 (?:(\d\d?(?:\.\d*)?)m)?\s*
887 (?:(\d\d?(?:\.\d*)?)s)?
888 \s*$ 'xi) {
889 $pord = $1 || "+";
890 $hd = $2 || 0;
891 # detect hord from format or from variable
892 $hord = ($3 && lc $3) || $hord;
893 $m = $4 || 0;
894 $s = $5 || 0;
895 $x = ($hd + $m/60. + $s/3600.) *
896 ($hord eq "h" ? 15 : 1) *
897 ($pord eq "-" ? -1 : 1);
898 } else { undef $x; } # Bad hms format
899 } else { # Not recognizeable
900 undef $x;
901 }
902
903 return $x;
904}
905
906# Convert a (lon,lat) pair in decimal degrees to sexagesimal representation.
907# Three decimal places of hours or two of degrees are kept.
908sub to_sex {
909 # v => hhmmss.ss instead of hhH mmM ss.sssS
910 my ($lon,$lat,$hord,$brief) = @_;
911
912 if($lon < 0) { $lon += 360.; }
913 (&dd_hdms($lon,3,$hord,$brief),&dd_hdms($lat,2,'d',$brief));
914}
915
916# Convert a single lon or lat in decimal hours or degrees to sexagesimal.
917# Number of decimal places on sec.s is adjustable.
918sub dd_hdms {
919 my($t,$dp,$hord,$brief) = @_;
920 my($h,$m,$s) = (0,0,0.0);
921 my($i,$sign,$str,$l,$cs);
922
923 $dp = defined $dp?$dp:3; $hord = lc ($hord || 'h');
924 $l = $dp+($dp?3:2);
925
926 $sign = $t =~ /\s*([+-])/ ? $1 : "" ;
927
928 $t = abs($t)/($hord eq 'h' ? 15 : 1)*3600.0;
929 $t = sprintf("%".
930 length(substr($t,0,(($i=index($t,"."))>=0?$i:length($t)))) .
931 ".${dp}f",$t);
932
933 $h = int $t/3600.0;
934 $m = int (($t - $h*3600.0)/60.);
935 $s = $t - ($h*3600.0 + $m*60.0);
936 $cs = sprintf("%0$l.${dp}f",$s);
937
938 if(! $brief) {
939 $str = sprintf("%s%02d%s%02dm${cs}s", $sign,$h,$hord,$m,$s);
940 } else {
941 $str = sprintf("%s%02d%02d${cs}", $sign,$h,$m,$s);
942 }
943
944 return $str;
945}
946
947# This is a copy of what's in Utils.pm, copied here to make
948# this module smaller if other modules ar enot otherwise needed
949
950# Make undefined values visible/printable
951sub _undefize {
952 my @x = map {defined $_ ? $_ : "<undef>"} @_;
953 return wantarray ? @x : $x[0];
954}
955
956# Regex to match a floating point number
957
# spent 37µs within WISE::CoUtils::_fpre which was called # once (37µs+0) at line 868
sub _fpre {
95872.7e-53.9e-6 my $d = shift || "";
959 $d = 'dD' if $d; # Is 'd' allowed instead of 'e' in the exponent?
960 my $re =
961 "(?:".
962 "[-+]?". # Optional_sign
963 "(?:".
964 "(?:(?:\\d*\\.?\\d+|\\d+\\.)(?:[eE$d][+-]?\\d+)?)|". # Normal_number
965 "(?:nanq?|inf(?:inity))". # Special_values
966 ")".
967 ")";
968 $re =~ s/nanq/[Nn][Aa][Nn][Qq]/; # Ensure case insenstivity w/o using (?i)
969 $re =~ s/inf/[Ii][Nn][Ff]/;
970 $re =~ s/inity/[Ii][Nn][Ii][Tt][Yy]/;
971 return $re;
972}
973
97413.0e-53.0e-51;

Report produced by the NYTProf 2.07 Perl profiler, developed by Tim Bunce and Adam Kaplan.