1package PDL::Graphics::Limits; 2 3use strict; 4use warnings; 5 6require Exporter; 7 8our @ISA = qw(Exporter); 9 10our %EXPORT_TAGS = ( 'all' => [ qw( 11 limits 12) ] ); 13 14our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } ); 15 16our @EXPORT = qw( 17 limits 18); 19 20our $VERSION = '0.01'; 21$VERSION = eval $VERSION; 22 23 24# Preloaded methods go here. 25 26use PDL::Core qw( cat pdl ); 27use PDL::Primitive qw( append ); 28use PDL::Fit::Polynomial; 29use PDL::Options; 30use PDL::Bad; 31use Carp; 32use POSIX qw( log10 ); 33 34use strict; 35use warnings; 36 37################################################################################ 38# figure out what's good in a piddle after a possible transformation which could 39# generate Infs or NaN's. If only everyone used PDL::Bad::UseNaN... 40sub set_mask 41{ 42 my ( $mask, $data ) = @_; 43 44 if ( $PDL::Bad::Status ) 45 { 46 my $badflag = $data->badflag(); 47 $data->badflag(1); 48 49 $mask .= $PDL::Bad::UseNaN ? (! $data->isbad ) : ( $data->isfinite & ! $data->isbad ); 50 51 $data->badflag($badflag); 52 } else 53 { 54 $mask .= $data->isfinite; 55 } 56} 57 58 59 60{ 61 package PDL::Graphics::Limits::DSet; 62 63 use PDL::Core qw( cat pdl ); 64 65 *set_mask = \*PDL::Graphics::Limits::set_mask; 66 67 sub new 68 { 69 my $class = shift; 70 my $self = bless {}, $class; 71 72 my ( $min, $max ) = splice( @_, 0, 2 ); 73 74 $self->{Vectors} = [ @_ ]; 75 $self->{MinMax} = [ map{ [ $min, $max ] } 1..@{$self->{Vectors}} ]; 76 77 $self; 78 } 79 80 sub ndim { scalar @{$_[0]->{Vectors}} } 81 82 sub validate 83 { 84 my ( $self, $attr) = @_; 85 86 my $ivec = 0; 87 my $n; 88 foreach my $vec ( @{$self->{Vectors}} ) 89 { 90 die( 'vector ', $ivec+1, ": no data?\n" ) 91 unless defined $vec->{data}; 92 93 $n = $vec->{data}->nelem unless defined $n; 94 95 # if a data set vector has no transformation function, use the 96 # default in $attr{Trans} 97 $vec->{trans} = $attr->{Trans}[$ivec] 98 if ! exists $vec->{trans} && exists $attr->{Trans}[$ivec]; 99 100 # remove explicitly undefined trans 101 delete $vec->{trans} 102 if exists $vec->{trans} && ! defined $vec->{trans}; 103 104 # ensure that data and errors have the same length. 105 die( 'vector ', $ivec+1, ": attribute $_: ", 106 "inconsistent number of elements", 107 "expected $n, got ", $vec->{$_}->nelem, "\n" ) 108 foreach 109 grep { exists $vec->{$_} && 110 defined $vec->{$_} && 111 $vec->{$_}->nelem != $n } 112 qw( data en ep ); 113 } 114 continue 115 { 116 $ivec++; 117 } 118 119 } 120 121 sub vector 122 { 123 $_[0]->{Vectors}[$_[1]]; 124 } 125 126 sub set_minmax 127 { 128 my ( $dset, $min, $max, $axis ) = @_; 129 130 my $mm = $dset->{MinMax}[$axis]; 131 132 $mm->[0] = $min if defined $min; 133 $mm->[1] = $max if defined $max; 134 } 135 136 sub upd_minmax 137 { 138 my ( $dset, $min, $max, $axis ) = @_; 139 140 my $mm = $dset->{MinMax}[$axis]; 141 142 $mm->[0] = $min if $mm->[0] > $min; 143 $mm->[1] = $max if $mm->[1] < $max; 144 } 145 146 sub get_minmax 147 { 148 my ( $dset ) = @_; 149 cat( map { pdl( $dset->{MinMax}[$_] ) } 0..$dset->ndim-1 ); 150 } 151 152 sub calc_minmax 153 { 154 my $dset = shift; 155 156 my @axes = @_ ? ( $_[0] ) : ( 0 ..$dset->ndims-1 ); 157 158 $dset->calc_minmax_axis( $_ ) foreach @axes; 159 } 160 161 ##################################################################### 162 # determine the limits for a dataset. 163 sub calc_minmax_axis 164 { 165 my ( $dset, $axis ) = @_; 166 167 my $vec = $dset->{Vectors}[$axis]; 168 my $data = $vec->{data}; 169 170 my $xfrm = defined $vec->{trans}; 171 172 # we need the transformed data point min max in case 173 # a transformed data + error is out of range of the transform 174 # function (e.g. log(0)). 175 176 my @minmax; 177 178 # reuse these as much as possible to reduce memory hit 179 my $tmp; 180 my $mask = PDL::null; 181 182 # i know of no way of determining whether a function can be applied inplace. 183 # assume not. 184 185 # if xfrm is true, $tmp will be an independent piddle, else its an alias for data 186 # no need to create a new piddle unless necessary. 187 $tmp = $xfrm ? $vec->{trans}->($data) : $data; 188 set_mask( $mask, $tmp ); 189 push @minmax, $tmp->where($mask)->minmax; 190 191 if ( defined $vec->{errn} ) 192 { 193 # worry about not overwriting the original data! 194 if ( $xfrm ) { $tmp .= $vec->{trans}->($data - $vec->{errn}) } 195 else { $tmp = $data - $vec->{errn} } 196 set_mask( $mask, $tmp ); 197 push @minmax, $tmp->where($mask)->minmax; 198 } 199 200 if ( defined $vec->{errp} ) 201 { 202 # worry about not overwriting the original data! 203 if ( $xfrm ) { $tmp .= $vec->{trans}->($data + $vec->{errp}) } 204 else { $tmp = $data + $vec->{errp} } 205 set_mask( $mask, $tmp ); 206 push @minmax, $tmp->where($mask)->minmax; 207 } 208 209 my ( $min, $max ) = PDL::Core::pdl( @minmax )->minmax; 210 211 $dset->set_minmax( $min, $max, $axis ); 212 } 213 214} 215 216##################################################################### 217 218# based upon PGPLOT's pgrnge routine. 219sub range_frac 220{ 221 my ( $axis, $frac, $zerofix ) = @_; 222 223 my $expand = $frac * ( $axis->[1] - $axis->[0] ); 224 my $min = $axis->[0] - $expand; 225 my $max = $axis->[1] + $expand; 226 227 if ( $zerofix ) 228 { 229 $min = 0.0 230 if $min < 0 && $axis->[0] >= 0.0; 231 232 $max = 0.0 233 if $max > 0 && $axis->[1] <= 0.0; 234 } 235 236 @{$axis} = ( $min, $max ); 237} 238 239##################################################################### 240 241# based upon PGPLOT's pgrnd routine 242 243# routine to find the closest "round" number to X, a "round" number 244# being 1, 2 or 5 times a power of 10. 245 246# If X is negative, round_pow(X) = -round_pow(abs(X)). 247# If X is zero, the value returned is zero. 248 249# round_pow( direction, $x ) 250# where direction is up, down, or both i.e. 251# $ub = round ( up => $x ); 252# $lb = round ( down => $x ); 253 254our @nice = ( 1, 2, 5, 10 ); 255our %flip = ( 'up' => 'down', 'down' => 'up' ); 256sub round_pow 257{ 258 my ( $what, $x ) = @_; 259 260 croak( "incorrect number of arguments" ) 261 unless 2 == @_; 262 263 if ( $x != 0.0 ) 264 { 265 my $xx = abs($x); 266 my $xlog = log10($xx); 267 my $ilog = int($xlog); 268 269 $what = $flip{$what} if $x < 0 ; 270 271 $ilog-- 272 if ( $xlog <= 0 && ( 'down' eq $what || $xlog != $ilog ) ) 273 || 274 ( $xlog > 0 && 'down' eq $what && $xlog == $ilog ) ; 275 276 my $pwr = 10 ** $ilog; 277 my $frac = $xx / $pwr; 278 279 my $i; 280 if ( 'up' eq $what ) 281 { 282 $i = 3; 283 $i = 2 if $frac < $nice[2]; 284 $i = 1 if $frac < $nice[1]; 285 $i = 0 if $frac < $nice[0]; 286 my $t = ( $x < 0 ? -1 : 1 ) * $pwr * $nice[$i]; 287 if(abs($t - $x) < 0.0000001) {$i++} 288 } 289 290 elsif ( 'down' eq $what ) 291 { 292 $i = 0; 293 $i = 1 if $frac > $nice[1]; 294 $i = 2 if $frac > $nice[2]; 295 $i = 3 if $frac > $nice[3]; 296 } 297 298 $x = ( $x < 0 ? -1 : 1 ) * $pwr * $nice[$i]; 299 } 300 301 $x; 302} 303 304##################################################################### 305 306sub setup_multi 307{ 308 my ( $common, $dim, $keys ) = @_; 309 310 my @arr; 311 312 if ( 'ARRAY' eq ref $common ) 313 { 314 return $common; 315 } 316 317 elsif ( 'HASH' eq ref $common ) 318 { 319 @arr[ 0..($dim-1)] = map { $common->{$_->{data}} } @{$keys}; 320 } 321 322 else 323 { 324 my $value = $common; 325 @arr = ($value) x $dim; 326 } 327 328 \@arr; 329} 330 331##################################################################### 332# normalize_dsets 333# 334# transform the user's heterogeneous list of data sets into a regular 335# list of data sets, each with the form 336# { Vectors => \@vectors } 337# where each vector is a hashref with the following keys: 338# { data => $data, 339# en => $err_n, 340# ep => $err_p, 341# trans => $trans } 342 343sub normalize_dsets 344{ 345 my ( $attr, @udsets ) = @_; 346 my @dsets; 347 348 while ( @udsets ) 349 { 350 my $ds = shift @udsets; 351 my $ref = ref $ds; 352 353 # peek inside the array to see what's there. we can have the following 354 # [ scalar|piddle, scalar|piddle, ... ] -> a zero dimensional data set 355 # [ \@a, \@b, \@c, \%d, ... ] -> a bunch of data sets 356 # [ \%h, @keys ] -> a hash with its keys 357 358 # scalar or piddle, turn it into its own data set 359 if ( ! $ref || UNIVERSAL::isa($ds, 'PDL') ) 360 { 361 push @dsets, 362 PDL::Graphics::Limits::DSet->new( $attr->{Min}, $attr->{Max}, 363 { data => PDL::Core::topdl( $ds ) } ); 364 } 365 366 elsif ( 'ARRAY' eq $ref ) 367 { 368 normalize_array( \@dsets, $attr, $ds ); 369 } 370 371 else 372 { 373 die( "data set: ", scalar @dsets + 1, 374 "illegal type in data set list: not an arrayref, scalar, or piddle\n" ); 375 } 376 377 } 378 379 # ensure data sets have the same dimensions 380 my %dim; 381 $dim{$_->ndim}++ foreach @dsets; 382 383 # whoops. only one allowed 384 die( "data sets do not all have the same dimensionality\n" ) 385 if keys %dim > 1; 386 387 ( $attr->{dims} ) = keys %dim; 388 389 # clean up datasets. 390 my $idset = -1; 391 foreach my $dset ( @dsets ) 392 { 393 $idset++; 394 395 eval { $dset->validate( $attr ) }; 396 if ( $@ ) 397 { 398 chomp $@; 399 die( "data set $idset: $@\n" ); 400 } 401 } 402 403 @dsets; 404} 405 406##################################################################### 407 408# array refs in data set lists may be just a plain ol' data set, or 409# it may contain a bunch of other stuff. here we deal with a single 410# array ref. we tear it apart and (re)build data sets. 411sub normalize_array 412{ 413 my ( $dsets, $attr, $aref ) = @_; 414 415 # if the first element is a hash, it's either a hash based data set 416 # with a bunch of attributes specific to that hash: 417 # [ \%h, @keys ] -> a hash with its keys 418 # in which case the rest of the elements are scalars, or its 419 # all hashes. 420 421 eval 422 { 423 if ( 'HASH' eq ref $aref->[0] ) 424 { 425 426 # all hashes? 427 if ( @$aref == grep { 'HASH' eq ref $_ } @$aref ) 428 { 429 # can't do anything unless we've been told which hash keys 430 # we should use, as this format doesn't allow local specification 431 die( "must specify hash keys for hash based data set spec\n" ) 432 unless defined $attr->{KeySpec} && scalar @{$attr->{KeySpec}}; 433 434 foreach ( @{$aref} ) 435 { 436 push @$dsets, normalize_hash_dset($attr, $_, @{$attr->{Keys}} ); 437 } 438 } 439 440 # hash + scalars? 441 elsif ( @$aref > 1 && 1 == grep { ref $_ } @$aref ) 442 { 443 push @$dsets, normalize_hash_dset( $attr, @{$aref} ) 444 } 445 446 # something wrong 447 else 448 { 449 die( "hash based data specification has an unexpected element" ); 450 } 451 452 } 453 454 # must be a list of vectors as either scalars, piddles, or array 455 # refs (vectors with attributes) 456 else 457 { 458 # for array based data sets, we have to accumulate vectors as we iterate 459 # through the array. they are stored here 460 my @vecs; 461 462 for my $vec ( @$aref ) 463 { 464 my $ref = ref $vec; 465 466 eval 467 { 468 # naked scalar or piddle: data vector with no attributes 469 if ( ! $ref || UNIVERSAL::isa($vec, 'PDL') ) 470 { 471 push @vecs, { data => PDL::Core::topdl( $vec ) }; 472 } 473 474 # array: data vector with attributes 475 elsif ( 'ARRAY' eq $ref ) 476 { 477 push @vecs, normalize_array_vec( $vec ); 478 } 479 480 else 481 { 482 die( 'vector ', @vecs+1, ": unexpected data type ($ref) in list of data sets\n" ); 483 } 484 }; 485 486 if ( $@ ) 487 { 488 chomp $@; 489 die( 'vector ', @vecs+1, ": $@\n" ); 490 } 491 } 492 493 push @$dsets, 494 PDL::Graphics::Limits::DSet->new( $attr->{Min}, $attr->{Max}, @vecs ) 495 if @vecs; 496 } 497 }; 498 499 if ( $@ ) 500 { 501 chomp $@; 502 die( 'data set ', @$dsets+1, ": $@\n" ); 503 } 504} 505 506##################################################################### 507 508# parse an array based vector 509sub normalize_array_vec 510{ 511 my ( $vec ) = @_; 512 513 # we should have 514 # [ $data, [ $err | $err_n, $err_p ], [ \&func ] ] 515 516 my @el = @$vec; 517 518 die( "too few or too many entries in array based data set spec\n" ) 519 if @el < 1 || @el > 4; 520 521 my %vec; 522 $vec{data} = PDL::Core::topdl( shift @el); 523 524 # if last value is CODE, it's a trans 525 $vec{trans} = pop @el if 'CODE' eq ref $el[-1]; 526 527 if ( exists $el[2] ) 528 { 529 # if we have 3 elements and the last isn't undef, it's an error. 530 # it can't be CODE as we'd have stripped it off in the last statement 531 die( "illegal value for trans func: $el[2]\n" ) 532 if defined $el[2]; 533 534 # we need to turn off trans for this element 535 $vec{trans} = undef; 536 pop @el; 537 } 538 539 # two values? asymmetric errors 540 if ( @el == 2 ) 541 { 542 $vec{errn} = PDL::Core::topdl($el[0]) if defined $el[0]; 543 $vec{errp} = PDL::Core::topdl($el[1]) if defined $el[1]; 544 } 545 546 # one value? symmetric errors 547 elsif ( @el == 1 ) 548 { 549 $vec{errn} = PDL::Core::topdl($el[0]) if defined $el[0]; 550 $vec{errp} = $vec{errn} if defined $vec{errn}; 551 } 552 553 \%vec; 554} 555 556##################################################################### 557 558# this takes a hash and a hash key spec and generates a regularized 559# data set array of the form 560# [ { data => $data, ep => ..., en => ..., trans => }, ... ] 561sub normalize_hash_dset 562{ 563 my ( $attr, $ds, @keys ) = @_; 564 565 my $KeySpec = $attr->{KeySpec}; 566 567 my @dset; 568 569 die( "too many local VecKeys (", scalar @keys, 570 ") and global VecKeys (", scalar @{$KeySpec}, ")\n" ) 571 if @keys && @{$KeySpec} && @{$KeySpec} <= @keys; 572 573 my @spec; 574 575 # handle local keys 576 if ( @keys ) 577 { 578 my $nvec = 0; 579 for my $key ( @keys ) 580 { 581 my %spec; 582 583 584 # parse the specs for this vector 585 eval { %spec = parse_vecspec( $key ) }; 586 do { chomp $@; die( "vector $nvec: $@" ) } 587 if $@; 588 589 590 # now, merge it with the global KeySpecs 591 592 if ( @{$KeySpec} ) 593 { 594 my $Spec = $KeySpec->[$nvec]; 595 596 foreach ( keys %{$Spec} ) 597 { 598 # only copy from Spec if not present in spec 599 $spec{$_} = $Spec->{$_} if ! exists $spec{$_}; 600 } 601 } 602 603 push @spec, \%spec; 604 } 605 continue 606 { 607 $nvec++; 608 } 609 610 # handle case where local VecKeys are a subst of global VecKeys 611 while ( @{$KeySpec} > @spec ) 612 { 613 push @spec, $KeySpec->[$nvec++]; 614 } 615 } 616 617 # no local keys; use global KeySpec 618 else 619 { 620 @spec = @{$KeySpec}; 621 } 622 623 my $nvec = 0; 624 for my $spec ( @spec ) 625 { 626 $nvec++; 627 my %vec; 628 629 die( "vector $nvec: no data spec?\n" ) 630 unless exists $spec->{data}; 631 632 for my $el ( qw( data errn errp trans ) ) 633 { 634 if ( exists $spec->{$el} ) 635 { 636 637 # if not defined, don't bother looking for it in the data set 638 unless ( defined $spec->{$el} ) 639 { 640 # trans is different from the others in that we need to pass 641 # it as undef if $spec->{trans} is undef (as full handling of 642 # trans is done elsewhere. 643 $vec{trans} = undef if 'trans' eq $el; 644 } 645 646 elsif ( exists $ds->{$spec->{$el}} ) 647 { 648 $vec{$el} = $ds->{$spec->{$el}}; 649 } 650 651 elsif ( $attr->{KeyCroak} ) 652 { 653 die( "vector $nvec: missing key in data set hash: ", $spec->{$el}, "\n" ) 654 } 655 } 656 657 } 658 659 # missing data; certainly a fatal error. 660 die( "vector $nvec: no data for key $spec->{data}\n" ) 661 unless defined $vec{data}; 662 663 push @dset, \%vec; 664 } 665 666 PDL::Graphics::Limits::DSet->new( $attr->{Min}, $attr->{Max}, @dset ); 667} 668 669##################################################################### 670# parse specifications for a hash based data set. These are the elements 671# in the VecKeys attribute. See the docs for more details. 672# Returns a hashref with keys data, en, ep, trans 673 674my $colre = qr/[^&<>=]/; 675 676# these are the different specs available. 677my %keyre = ( data => qr/^($colre+)/, 678 errn => qr/<($colre*)/, 679 errp => qr/>($colre*)/, 680 err => qr/=($colre*)/, 681 trans => qr/\&($colre*)/ 682 ); 683 684my %vecspeckeys = ( data => 1, 685 err => 1, 686 errn => 1, 687 errp => 1, 688 trans => 1 ); 689 690sub parse_vecspec 691{ 692 my ( $ukeys ) = @_; 693 694 my %k; 695 696 # do we get a hash? 697 if ( 'HASH' eq ref $ukeys ) 698 { 699 # complain about keys we don't use 700 my @badkeys = grep { ! defined $vecspeckeys{$_} } keys %$ukeys; 701 die( "illegal keys: ", join(' ,', @badkeys), "\n" ) 702 if @badkeys; 703 704 # copy keys we need 705 do { $k{$_} = $ukeys->{$_} if exists $ukeys->{$_} } 706 foreach keys %vecspeckeys; 707 708 } 709 710 # parse the string. 711 else 712 { 713 714 # make a local copy, as we modify it in place. 715 my $keys = $ukeys; 716 717 # deal with a "default" spec 718 if ( ! defined $keys ) 719 { 720 $keys = ''; 721 } 722 else 723 { 724 # spaces and commas are there for human use only 725 $keys =~ s/[\s,]//g; 726 } 727 728 729 # extract the known specs. 730 my ( $what, $re ); 731 $keys =~ s/$re// and $k{$what} = $1 while( ($what, $re) = each %keyre); 732 733 # if there's anything left, it's bogus 734 die( "illegal key specification: $ukeys\n" ) 735 unless $keys eq ''; 736 737 } 738 739 # check for consistent error bar specs 740 die( "can't specify `=' with `<' or `>'\n" ) 741 if exists $k{err} && ( exists $k{errn} || exists $k{errp} ); 742 743 # error bars are always specified as positive and negative; turn a symmetric 744 # spec into that 745 $k{errn} = $k{errp} = $k{err} if exists $k{err}; 746 delete $k{err}; 747 748 # set empty values to undefined ones 749 do { $k{$_} = undef if $k{$_} eq '' } foreach keys %k; 750 751 %k; 752} 753 754sub parse_vecspecs 755{ 756 my $keys = shift; 757 my @specs; 758 759 push @specs, { parse_vecspec($_) } 760 foreach @$keys; 761 762 \@specs; 763} 764 765##################################################################### 766# normalize user supplied limits 767 768sub parse_limits 769{ 770 my ( $ndim, $spec, $KeySpec ) = @_; 771 772 $spec = [] unless defined $spec; 773 774 my @limits; 775 776 # array containing limits (as arrays or scalars) 777 if ( 'ARRAY' eq ref $spec ) 778 { 779 # no limits; just move on 780 unless ( @$spec ) 781 { 782 } 783 784 # multi-dimensional data sets 785 elsif ( 'ARRAY' eq ref $spec->[0] ) 786 { 787 my $ilim = 0; 788 for my $vlim ( @$spec ) 789 { 790 $ilim++; 791 die( "Limit spec element $ilim: expected array ref\n" ) 792 if 'ARRAY' ne ref $vlim; 793 794 die( "Limit spec element $ilim: too many values\n" ) 795 if @$vlim > 2; 796 797 die( "Limit spec element $vlim: values must be scalars\n" ) 798 if grep { ref $_ } @$vlim; 799 800 my @lims = @$vlim; 801 $lims[0] = undef unless defined $lims[0]; 802 $lims[1] = undef unless defined $lims[1]; 803 804 push @limits, \@lims; 805 } 806 } 807 808 # one-dimensional data sets 809 elsif ( ! ref $spec->[0] ) 810 { 811 die( "unexpected non-scalar element in Limits spec\n" ) 812 if grep { ref $_ } @$spec; 813 814 my @lims = @$spec; 815 $lims[0] = undef unless defined $lims[0]; 816 $lims[1] = undef unless defined $lims[1]; 817 818 push @limits, \@lims; 819 } 820 821 push @limits, [ undef, undef ] 822 while ( @limits != $ndim ); 823 824 } 825 826 # hash containing vector names and limits 827 elsif ( 'HASH' eq ref $spec ) 828 { 829 # first ensure that VecKeys has been specified 830 die( "cannot use Limits without VecKeys\n" ) 831 unless @$KeySpec; 832 833 # make sure that we've got common keys. 834 my %vecs = map { ( $_->{data} => 1) } @$KeySpec; 835 836 # identify unknown vectors 837 my @badvecs = grep { ! defined $vecs{$_} } keys %$spec; 838 die( 'unknown vector(s): ', join(', ', @badvecs), "\n" ) 839 if @badvecs; 840 841 # work our way through the KeySpec's, filling in values from 842 # $spec as appropriate. 843 for my $kspec ( @$KeySpec ) 844 { 845 my @lims = ( undef, undef ); 846 if ( exists $spec->{$kspec->{data}} ) 847 { 848 my $lspec = $spec->{$kspec->{data}}; 849 $lims[0] = $lspec->{min} if exists $lspec->{min}; 850 $lims[1] = $lspec->{max} if exists $lspec->{max}; 851 } 852 push @limits, \@lims; 853 } 854 } 855 856 # say what? 857 else 858 { 859 die( "Limits attribute value must be a hashref or arrayref\n" ); 860 } 861 862 map { { calc => scalar ( grep { !defined $_ } @{$_} ), range => $_ } } @limits; 863} 864 865 866 867##################################################################### 868 869sub limits 870{ 871 my $attr = 'HASH' eq ref $_[-1] ? pop @_ : {}; 872 873 my @udsets = @_; 874 875 my %attr = iparse( { 876 Min => -1.8e308, 877 Max => +1.8e308, 878 Bounds => 'minmax', 879 Clean => 'RangeFrac', 880 RangeFrac => 0.05, 881 ZeroFix => 0, 882 VecKeys => [], 883 KeyCroak => 1, 884 Limits => [], 885 Trans => [], 886 }, $attr ); 887 888 # turn Trans and VecKeys into arrays if necessary; may be scalars for 1D 889 # data sets 890 $attr{$_} = [ $attr{$_} ] 891 foreach grep { defined $attr{$_} && 'ARRAY' ne ref $attr{$_} } 892 qw( VecKeys Trans ); 893 894 # parse vector key specs 895 $attr{KeySpec} = parse_vecspecs( $attr{VecKeys} ); 896 897 # normalize data sets to make life easier later. also 898 # counts up the number of dimensions and sets $attr{dims} 899 my @dsets = normalize_dsets( \%attr, @udsets ); 900 901 # set up the Limits 902 my @limits = parse_limits( $attr{dims}, $attr{Limits}, $attr{KeySpec} ); 903 904 if ( 'minmax' eq lc $attr{Bounds} ) 905 { 906 for my $dim ( 0..$attr{dims}-1 ) 907 { 908 # only calculate minmax values for those dims which need them. 909 my $limits = $limits[$dim]; 910 911 foreach ( @dsets ) 912 { 913 # calculate min & max 914 $_->calc_minmax( $dim ) 915 if $limits->{calc}; 916 917 # make sure we pay attention to user specified limits 918 $_->set_minmax( @{$limits->{range}}, $dim ); 919 } 920 } 921 } 922 923 elsif ( 'zscale' eq lc $attr{Bounds} ) 924 { 925 croak( "zscale only good for dim = 2\n" ) 926 unless $attr{dims} == 2; 927 928 foreach my $dset ( @dsets ) 929 { 930 $dset->calc_minmax( 0 ) 931 if $limits[0]{calc}; 932 933 934 if ( $limits[1]{calc} ) 935 { 936 my $y = $dset->vector(1)->{data}; 937 938 # this is a waste, as we don't care about the evaluated 939 # fit values, just the min and max values. since we 940 # get them all anyway, we'll use them. 941 942 my $mask = PDL::null; 943 set_mask( $mask, $y ); 944 945 my $fit = fitpoly1d( $y->where($mask)->qsort, 2 ); 946 $dset->set_minmax( $fit->minmax, 1 ); 947 } 948 949 $dset->set_minmax( @{$limits[$_]{range}}, $_ ) for 0,1; 950 } 951 } 952 else 953 { 954 die( "unknown Bounds type: $attr{Bounds}\n" ); 955 } 956 957 # derive union of minmax limits from data sets 958 my $minmax = PDL::Core::null; 959 $minmax = append( $minmax, $_->get_minmax ) foreach @dsets; 960 961 # get overall minmax limits 962 $minmax = cat(($minmax->minmaximum)[0,1])->transpose; 963 964 my @minmax = map{ [ $minmax->slice(":,$_")->list ] } 0..$attr{dims}-1; 965 966 if ( 'rangefrac' eq lc $attr{Clean} ) 967 { 968 my $RangeFrac = 969 setup_multi( $attr{RangeFrac}, $attr{dims}, $attr{KeySpec} ); 970 971 my $ZeroFix = 972 setup_multi( $attr{ZeroFix}, $attr{dims}, $attr{KeySpec} ); 973 974 range_frac( $minmax[$_], $RangeFrac->[$_], $ZeroFix->[$_] ) 975 for 0..$attr{dims}-1; 976 } 977 978 elsif ( 'roundpow' eq lc $attr{Clean} ) 979 { 980 $_ = [ round_pow( down => $_->[0] ), 981 round_pow( up => $_->[1] ) ] 982 foreach @minmax; 983 } 984 985 elsif ( 'none' eq lc $attr{Clean} ) 986 { 987 # do nothing 988 } 989 990 else 991 { 992 die( "unknown Clean type: $attr{Clean}\n" ); 993 } 994 995 if ( wantarray ) 996 { 997 return map { ( @{$_} ) } @minmax; 998 } 999 else 1000 { 1001 my @key; 1002 if ( @{$attr{KeySpec}} ) 1003 { 1004 @key = map { $_->{data} } @{$attr{KeySpec}}; 1005 } 1006 else 1007 { 1008 @key = map { 'q' . $_ } ( 1 .. $attr{dims} ); 1009 } 1010 1011 return { map { ( $key[$_] => { min => $minmax[$_][0], 1012 max => $minmax[$_][1] } ) } 1013 0.. ( @minmax - 1 ) }; 1014 } 1015} 1016 10171; 1018 1019 1020__END__ 1021 1022=pod 1023 1024=head1 NAME 1025 1026PDL::Graphics::Limits - derive limits for display purposes 1027 1028 1029=head1 DESCRIPTION 1030 1031Functions to derive limits for data for display purposes 1032 1033=head1 SYNOPSIS 1034 1035 use PDL::Graphics::Limits; 1036 1037 1038=head1 FUNCTIONS 1039 1040=head2 limits 1041 1042=for ref 1043 1044B<limits> derives global limits for one or more multi-dimensional sets 1045of data for display purposes. It obtains minimum and maximum limits 1046for each dimension based upon one of several algorithms. 1047 1048=for usage 1049 1050 @limits = limits( @datasets ); 1051 @limits = limits( @datasets, \%attr ); 1052 $limits = limits( @datasets ); 1053 $limits = limits( @datasets, \%attr ); 1054 1055=head3 Data Sets 1056 1057A data set is represented as a set of one dimensional vectors, one per 1058dimension. All data sets must have the same dimensions. 1059Multi-dimensional data sets are packaged as arrays or hashs; one 1060dimensional data sets need not be. The different representations may 1061be mixed, as long as the dimensions are presented in the same order. 1062Vectors may be either scalars or piddles. 1063 1064=over 8 1065 1066=item One dimensional data sets 1067 1068One dimensional data sets may be passed directly, with no additional packaging: 1069 1070 limits( $scalar, $piddle ); 1071 1072=item Data sets as arrays 1073 1074If the data sets are represented by arrays, each vectors in each array 1075must have the same order: 1076 1077 @ds1 = ( $x1_pdl, $y1_pdl ); 1078 @ds2 = ( $x2_pdl, $y2_pdl ); 1079 1080They are passed by reference: 1081 1082 limits( \@ds1, \@ds2 ); 1083 1084=item Data sets as hashes 1085 1086Hashes are passed by reference as well, but I<must> be further 1087embedded in arrays (also passed by reference), in order that the last 1088one is not confused with the optional trailing attribute hash. For 1089example: 1090 1091 limits( [ \%ds4, \%ds5 ], \%attr ); 1092 1093If each hash uses the same keys to identify the data, the keys 1094should be passed as an ordered array via the C<VecKeys> attribute: 1095 1096 limits( [ \%h1, \%h2 ], { VecKeys => [ 'x', 'y' ] } ); 1097 1098If the hashes use different keys, each hash must be accompanied by an 1099ordered listing of the keys, embedded in their own anonymous array: 1100 1101 [ \%h1 => ( 'x', 'y' ) ], [ \%h2 => ( 'u', 'v' ) ] 1102 1103Keys which are not explicitly identified are ignored. 1104 1105=back 1106 1107=head3 Errors 1108 1109Error bars must be taken into account when determining limits; care 1110is especially needed if the data are to be transformed before plotting 1111(for logarithmic plots, for example). Errors may be symmetric (a single 1112value indicates the negative and positive going errors for a data point) or 1113asymmetric (two values are required to specify the errors). 1114 1115If the data set is specified as an array of vectors, vectors with 1116errors should be embedded in an array. For symmetric errors, the error 1117is given as a single vector (piddle or scalar); for asymmetric errors, there 1118should be two values (one of which may be C<undef> to indicate 1119a one-sided error bar): 1120 1121 @ds1 = ( $x, # no errors 1122 [ $y, $yerr ], # symmetric errors 1123 [ $z, $zn, $zp ], # asymmetric errors 1124 [ $u, undef, $up ], # one-sided error bar 1125 [ $v, $vn, undef ], # one-sided error bar 1126 ); 1127 1128If the data set is specified as a hash of vectors, the names of the 1129error bar keys are appended to the names of the data keys in the 1130C<VecKeys> designations. The error bar key names are always prefixed 1131with a character indicating what kind of error they represent: 1132 1133 < negative going errors 1134 > positive going errors 1135 = symmetric errors 1136 1137(Column names may be separated by commas or white space.) 1138 1139For example, 1140 1141 %ds1 = ( x => $x, xerr => $xerr, y => $y, yerr => $yerr ); 1142 limits( [ \%ds1 ], { VecKeys => [ 'x =xerr', 'y =yerr' ] } ); 1143 1144To specify asymmetric errors, specify both the negative and positive going 1145errors: 1146 1147 %ds1 = ( x => $x, xnerr => $xn, xperr => $xp, 1148 y => $y ); 1149 limits( [ \%ds1 ], { VecKeys => [ 'x <xnerr >xperr', 'y' ] } ); 1150 1151For one-sided error bars, specify a column just for the side to 1152be plotted: 1153 1154 %ds1 = ( x => $x, xnerr => $xn, 1155 y => $y, yperr => $yp ); 1156 limits( [ \%ds1 ], { VecKeys => [ 'x <xnerr', 'y >yperr' ] } ); 1157 1158 1159Data in hashes with different keys follow the same paradigm: 1160 1161 [ \%h1 => ( 'x =xerr', 'y =yerr' ) ], [ \%h2 => ( 'u =uerr', 'v =verr' ) ] 1162 1163In this case, the column names specific to a single data set override 1164those specified via the C<VecKeys> option. 1165 1166 limits( [ \%h1 => 'x =xerr' ], { VecKeys => [ 'x <xn >xp' ] } ) 1167 1168In the case of a multi-dimensional data set, one must specify 1169all of the keys: 1170 1171 limits( [ \%h1 => ( 'x =xerr', 'y =yerr' ) ], 1172 { VecKeys => [ 'x <xn >xp', 'y <yp >yp' ] } ) 1173 1174One can override only parts of the specifications: 1175 1176 limits( [ \%h1 => ( '=xerr', '=yerr' ) ], 1177 { VecKeys => [ 'x <xn >xp', 'y <yp >yp' ] } ) 1178 1179Use C<undef> as a placeholder for those keys for which 1180nothing need by overridden: 1181 1182 limits( [ \%h1 => undef, 'y =yerr' ], 1183 { VecKeys => [ 'x <xn >xp', 'y <yp >yp' ] } ) 1184 1185=head3 Data Transformation 1186 1187Normally the data passed to B<limits> should be in their final, 1188transformed, form. For example, if the data will be displayed on a 1189logarithmic scale, the logarithm of the data should be passed to 1190B<limits>. However, if error bars are also to be displayed, the 1191I<untransformed> data must be passed, as 1192 1193 log(data) + log(error) != log(data + error) 1194 1195Since the ranges must be calculated for the transformed values, 1196B<range> must be given the transformation function. 1197 1198If all of the data sets will undergo the same transformation, this may 1199be done with the B<Trans> attribute, which is given a list of 1200subroutine references, one for each element of a data set. An 1201C<undef> value may be used to indicate no transformation is to be 1202performed. For example, 1203 1204 @ds1 = ( $x, $y ); 1205 1206 # take log of $x 1207 limits( \@ds1, { trans => [ \&log10 ] } ); 1208 1209 # take log of $y 1210 limits( \@ds1, { trans => [ undef, \&log10 ] } ); 1211 1212If each data set has a different transformation, things are a bit more 1213complicated. If the data sets are specified as arrays of vectors, vectors 1214with transformations should be embedded in an array, with the I<last> 1215element the subroutine reference: 1216 1217 @ds1 = ( [ $x, \&log10 ], $y ); 1218 1219With error bars, this looks like this: 1220 1221 @ds1 = ( [ $x, $xerr, \&log10 ], $y ); 1222 @ds1 = ( [ $x, $xn, $xp, \&log10 ], $y ); 1223 1224If the C<Trans> attribute is used in conjunction with individual data 1225set transformations, the latter will override it. To explicitly 1226indicate that a specific data set element has no transformation 1227(normally only needed if C<Trans> is used to specify a default) set 1228the transformation subroutine reference to C<undef>. In this case, 1229the entire quad of data element, negative error, positive error, and 1230transformation subroutine must be specified to avoid confusion: 1231 1232 [ $x, $xn, $xp, undef ] 1233 1234Note that $xn and $xp may be undef. For symmetric errors, simply 1235set both C<$xn> and C<$xp> to the same value. 1236 1237For data sets passed as hashes, the subroutine reference is an element 1238in the hashes; the name of the corresponding key is added to the list 1239of keys, preceded by the C<&> character: 1240 1241 %ds1 = ( x => $x, xerr => $xerr, xtrans => \&log10, 1242 y => $y, yerr => $yerr ); 1243 1244 limits( [ \%ds1, \%ds2 ], 1245 { VecKeys => [ 'x =xerr &xtrans', 'y =yerr' ] }); 1246 limits( [ \%ds1 => 'x =xerr &xtrans', 'y =yerr' ] ); 1247 1248If the C<Trans> attribute is specified, and a key name is also 1249specified via the C<VecKeys> attribute or individually for a data set 1250element, the latter will take precedence. For example, 1251 1252 $ds1{trans1} = \&log10; 1253 $ds1{trans2} = \&sqrt; 1254 1255 # resolves to exp 1256 limits( [ \%ds1 ], { Trans => [ \&exp ] }); 1257 1258 # resolves to sqrt 1259 limits( [ \%ds1 ], { Trans => [ \&exp ], 1260 VecKeys => [ 'x =xerr &trans2' ] }); 1261 1262 # resolves to log10 1263 limits( [ \%ds1 => '&trans1' ], { Trans => [ \&exp ], 1264 VecKeys => [ 'x =xerr &trans2' ] }); 1265 1266 1267To indicate that a particular vector should have no transformation, 1268use a blank key: 1269 1270 limits( [ \%ds1 => ( 'x =xerr &', 'y =yerr' ) ], [\%ds2], 1271 { Trans => [ \&log10 ] } ); 1272 1273or set the hash element to C<undef>: 1274 1275 $ds1{xtrans} = undef; 1276 1277 1278=head3 Range Algorithms 1279 1280Sometimes all you want is to find the minimum and maximum values. However, 1281for display purposes, it's often nice to have "clean" range bounds. To that 1282end, B<limits> produces a range in two steps. First it determines the bounds, 1283then it cleans them up. 1284 1285To specify the bounding algorithm, set the value of the C<Bounds> key 1286in the C<%attr> hash to one of the following values: 1287 1288=over 8 1289 1290=item MinMax 1291 1292This indicates the raw minima and maxima should be used. This is the 1293default. 1294 1295=item Zscale 1296 1297This is valid for two dimensional data only. The C<Y> values are sorted, 1298then fit to a line. The minimum and maximum values of the evaluated 1299line are used for the C<Y> bounds; the raw minimum and maximum values 1300of the C<X> data are used for the C<X> bounds. This method is good 1301in situations where there are "spurious" spikes in the C<Y> data which 1302would generate too large a dynamic range in the bounds. (Note that 1303the C<Zscale> algorithm is found in IRAF and DS9; its true origin 1304is unknown to the author). 1305 1306=back 1307 1308To specify the cleaning algorithm, set the value of the C<Clean> key 1309in the C<%attr> hash to one of the following values: 1310 1311=over 8 1312 1313=item None 1314 1315Perform no cleaning of the bounds. 1316 1317=item RangeFrac 1318 1319This is based upon the C<PGPLOT> B<pgrnge> function. It symmetrically expands 1320the bounds (determined above) by a fractional amount: 1321 1322 $expand = $frac * ( $axis->{max} - $axis->{min} ); 1323 $min = $axis->{min} - $expand; 1324 $max = $axis->{max} + $expand; 1325 1326The fraction may be specified in the C<%attr> hash with the 1327C<RangeFrac> key. It defaults to C<0.05>. 1328 1329Because this is a symmetric expansion, a limit of C<0.0> may be 1330transformed into a negative number, which may be inappropriate. If 1331the C<ZeroFix> key is set to a non-zero value in the C<%attr> hash, 1332the cleaned boundary is set to C<0.0> if it is on the other side of 1333C<0.0> from the above determined bounds. For example, If the minimum 1334boundary value is C<0.1>, and the cleaned boundary value is C<-0.1>, 1335the cleaned value will be set to C<0.0>. Similarly, if the maximum 1336value is C<-0.1> and the cleaned value is C<0.1>, it will be set to C<0.0>. 1337 1338This is the default clean algorithm. 1339 1340=item RoundPow 1341 1342This is based upon the C<PGPLOT> B<pgrnd> routine. It determines a 1343"nice" value, where "nice" is the closest round number to 1344the boundary value, where a round number is 1, 2, or 5 times a power 1345of 10. 1346 1347=back 1348 1349=head3 User Specified Limits 1350 1351To fully or partially override the automatically determined limits, 1352use the B<Limits> attribute. These values are used as input to the 1353range algorithms. 1354 1355The B<Limits> attribute value may be either an array of arrayrefs, or 1356a hash. 1357 1358=over 1359 1360=item Arrays 1361 1362The B<Limits> value may be a reference to an array of arrayrefs, one 1363per dimension, which contain the requested limits. 1364 1365The dimensions should be ordered in the same way as the datasets. 1366Each arrayref should contain two ordered values, the minimum and 1367maximum limits for that dimension. The limits may have the undefined 1368value if that limit is to be automatically determined. The limits 1369should be transformed (or not) in the same fashion as the data. 1370 1371For example, to specify that the second dimension's maximum limit 1372should be fixed at a specified value: 1373 1374 Limits => [ [ undef, undef ], [ undef, $max ] ] 1375 1376Note that placeholder values are required for leading dimensions which 1377are to be handled automatically. For convenience, if limits for a 1378dimension are to be fully automatically determined, the placeholder 1379arrayref may be empty. Also, trailing undefined limits may be 1380omitted. The above example may be rewritten as: 1381 1382 Limits => [ [], [ undef, $max ] ] 1383 1384If the minimum value was specified instead of the maximum, the following 1385would be acceptable: 1386 1387 Limits => [ [], [ $min ] ] 1388 1389If the data has but a single dimension, nested arrayrefs are not required: 1390 1391 Limits => [ $min, $max ] 1392 1393 1394=item Hashes 1395 1396Th B<Limits> attribute value may be a hash; this can only be used in 1397conjunction with the B<VecKeys> attribute. If the data sets are 1398represented by hashes which do not have common keys, then the user 1399defined limits should be specified with arrays. The keys in the 1400B<Limits> hash should be the names of the data vectors in the 1401B<VecKeys>. Their values should be hashes with keys C<min> and C<max>, 1402representing the minimum and maximum limits. Limits which have the value 1403C<undef> or which are not specified will be determined from the data. 1404For example, 1405 1406 Limits => { x => { min => 30 }, y => { max => 22 } } 1407 1408=back 1409 1410=head3 Return Values 1411 1412When called in a list context, it returns the minimum and maximum 1413bounds for each axis: 1414 1415 @limits = ( $min_1, $max_1, $min_2, $max_2, ... ); 1416 1417which makes life easier when using the B<env> method: 1418 1419 $window->env( @limits ); 1420 1421When called in a scalar context, it returns a hashref with the keys 1422 1423 axis1, ... axisN 1424 1425where C<axisN> is the name of the Nth axis. If axis names have not 1426been specified via the C<VecKeys> element of C<%attr>, names are 1427concocted as C<q1>, C<q2>, etc. The values are hashes with keys 1428C<min> and C<max>. For example: 1429 1430 { q1 => { min => 1, max => 2}, 1431 q2 => { min => -33, max => 33 } } 1432 1433=head3 Miscellaneous 1434 1435Normally B<limits> complains if hash data sets don't contain specific 1436keys for error bars or transformation functions. If, however, 1437you'd like to specify default values using the C<%attr> argument, 1438but there are data sets which don't have the data and you'd rather 1439not have to explicitly indicate that, set the C<KeyCroak> attribute 1440to zero. For example, 1441 1442 limits( [ { x => $x }, { x => $x1, xerr => $xerr } ], 1443 { VecKeys => [ 'x =xerr' ] } ); 1444 1445will generate an error because the first data set does not have 1446an C<xerr> key. Resetting C<KeyCroak> will fix this: 1447 1448 limits( [ { x => $x }, { x => $x1, xerr => $xerr } ], 1449 { VecKeys => [ 'x =xerr' ], KeyCroak => 0 } ); 1450 1451 1452=head1 AUTHOR 1453 1454Diab Jerius, E<lt>djerius@cpan.orgE<gt> 1455 1456=head1 COPYRIGHT AND LICENSE 1457 1458Copyright (C) 2004 by the Smithsonian Astrophysical Observatory 1459 1460 1461This software is released under the GNU General Public License. 1462You may find a copy at L<http://www.fsf.org/copyleft/gpl.html>. 1463 1464=cut 1465