NAME
PDL::Ufunc - primitive ufunc operations for pdl
DESCRIPTION
This module provides some primitive and useful functions defined using PDL::PP based on functionality of what are sometimes called ufuncs (for example NumPY and Mathematica talk about these). It collects all the functions generally used to reduce or accumulate along a dimension. These all do their job across the first dimension but by using the slicing functions you can do it on any dimension.
The PDL::Reduce module provides an alternative interface to many of the functions in this module.
SYNOPSIS
use PDL::Ufunc;Project via $name to N-1 dimensions
This function reduces the dimensionality of an ndarray by one by taking the $name along the 1st dimension.
By using xchg etc. it is possible to use any dimension.
\$y = $op(\$x); \$spectrum = $op \$image->transpose$extras
Cumulative $name
This function calculates the cumulative $name along the 1st dimension.
By using xchg etc. it is possible to use any dimension.
The sum is started so that the first element in the cumulative $name is the first element of the parameter.
\$y = $op(\$x); \$spectrum = $op \$image->transpose$extras
Numerically differentiates a vector along 0th dimension.
By using "xchg" in PDL::Slices etc. it is possible to use any dimension.
print pdl(q[3 4 2 3 2 3 1])->diff2;
# [1 -2 1 -1 1 -2]
EOF
BadDoc => <<'EOF',
On bad value, output value is set bad. On next good value, output value
is difference between that and last good value.
EOF
);# this would need a lot of work to support bad values # plus it gives me a chance to check out HandleBad => 0 ;) # pp_def( 'intover', HandleBad => 0, Pars => 'a(n); float+ [o]b();', Code => '/* Integration formulae from Press et al 2nd Ed S 4.1 */ switch ($SIZE(n)) { case 1: broadcastloop %{ $b() = 0.; /* not a(n=>0); as interval has zero width */ %} break; case 2: broadcastloop %{ $b() = 0.5*($a(n=>0)+$a(n=>1)); %} break; case 3: broadcastloop %{ $b() = ($a(n=>0)+4*$a(n=>1)+$a(n=>2))/3.; %} break; case 4: broadcastloop %{ $b() = ($a(n=>0)+$a(n=>3)+3.*($a(n=>1)+$a(n=>2)))*0.375; %} break; case 5: broadcastloop %{ $b() = (14.*($a(n=>0)+$a(n=>4)) +64.*($a(n=>1)+$a(n=>3)) +24.*$a(n=>2))/45.; %} break; default: broadcastloop %{ $GENERIC(b) tmp = 0; loop (n=3:-3) %{ tmp += $a(); %} loop (n=-3:-2) %{ tmp += (23./24.)*($a(n=>2)+$a()); %} loop (n=-2:-1) %{ tmp += (7./6.) *($a(n=>1)+$a()); %} loop (n=-1:) %{ tmp += (3./8.) *($a(n=>0)+$a()); %} $b() = tmp; %} } ', Doc => projectdocs('integral','intover', <<'EOF'), Notes:
intover uses a point spacing of one (i.e., delta-h==1). You will need to scale the result to correct for the true point delta.
For n > 3, these are all O(h^4) (like Simpson's rule), but are integrals between the end points assuming the pdl gives values just at these centres: for such `functions', sumover is correct to O(h), but is the natural (and correct) choice for binned data, of course. EOF ); # intover
sub synonym { my ($name, $synonym) = @_; pp_add_exported('', $synonym); pp_addpm( "=head2 $synonym\n\n=for ref\n\nSynonym for "$name".\n\n=cut\n *PDL::$synonym = *$synonym = \\&PDL::$name;" ); }
sub make_average { my ($prefix, $outpar_type, $extra) = @_; pp_def( "${prefix}average", HandleBad => 1, Pars => "a(n); $outpar_type [o]b();", Code => q{ $GENERIC(b) tmp = 0; PDL_Indx cnt = 0; loop(n) %{ PDL_IF_BAD(if ( $ISGOOD(a()) ) {,) cnt++; tmp += $a(); PDL_IF_BAD(},) %} if ( !cnt ) { PDL_IF_BAD($SETBAD(b()), $b() = PDL_IF_GENTYPE_INTEGER(0,NAN)); } else $b() = tmp / cnt; }, Doc => projectdocs( 'average', "${prefix}average", $extra||'' ), ); synonym(map "$prefix$_", qw(average avgover)); }
make_average('', 'int+'); make_average('c', 'cdouble', "Unlike average, the calculation is performed in complex double precision." ); make_average('d', 'double', "Unlike average, the calculation is performed in double precision." );
for my $which ( [qw(minimum < minover)], [qw(maximum > maxover)], ) { my ($name, $op, $synonym) = @$which;
pp_def(
$name,
HandleBad => 1,
Pars => 'a(n); [o]c();',
Code =>
'$GENERIC() cur = 0;
int flag = 0;
loop(n) %{
PDL_IF_BAD(if ($ISBAD(a())) continue;,)
if( !flag || ($a() '.$op.' cur ) || PDL_ISNAN_$PPSYM()(cur) )
{cur = $a(); flag = 1;}
%}
if ( flag ) { $c() = cur; }
else { $SETBAD(c()); $PDLSTATESETBAD(c); }',
Doc => projectdocs($name,$name,''),
BadDoc =>
'Output is set bad if no elements of the input are non-bad,
otherwise the bad flag is cleared for the output ndarray.Note that NaNs are considered to be valid values and will "win" over non-NaN; see isfinite and badmask for ways of masking NaNs. ', ); synonym($name, $synonym);
pp_def(
"${name}_ind",
HandleBad => 1,
Pars => 'a(n); indx [o] c();',
Code =>
'$GENERIC() cur = 0;
PDL_Indx curind = -1;
loop(n) %{
PDL_IF_BAD(if ($ISBAD(a())) continue;,)
if(curind == -1 || $a() '.$op.' cur || PDL_ISNAN_$PPSYM()(cur))
{cur = $a(); curind = n;}
%}
if ( curind != -1 ) { $c() = curind; }
else { $SETBAD(c()); $PDLSTATESETBAD(c); }',
Doc => "Like $name but returns the index rather than the value",
BadDoc =>
'Output is set bad if no elements of the input are non-bad,
otherwise the bad flag is cleared for the output ndarray.Note that NaNs are considered to be valid values and will "win" over non-NaN; see isfinite and badmask for ways of masking NaNs. ', ); synonym("${name}_ind", "${synonym}_ind");
pp_def(
"${name}_n_ind",
HandleBad => 1,
Pars => 'a(n); indx [o]c(m);',
OtherPars => 'PDL_Indx m_size => m;',
PMCode => PDL::PP::pp_line_numbers(__LINE__, <<EOF),
sub PDL::${name}_n_ind {
my (\$a, \$c, \$m_size) = \@_;
\$m_size //= ref(\$c) ? \$c->dim(0) : \$c; # back-compat with pre-2.077
my \$set_out = 1;
\$set_out = 0, \$c = null if !ref \$c;
\$c = \$c->indx if !\$c->isnull;
PDL::_${name}_n_ind_int(\$a, \$c, \$m_size);
\$set_out ? \$_[1] = \$c : \$c;
}
EOF
RedoDimsCode => 'if($SIZE(m) > $SIZE(n)) $CROAK("m_size > n_size");',
Code =>
'$GENERIC() cur = 0; PDL_Indx curind; register PDL_Indx ns = $SIZE(n);
$PDLSTATESETGOOD(c);
loop(m) %{
curind = ns;
loop(n) %{
PDL_Indx outer_m = m; int flag=0;
loop (m=:outer_m) %{
if ($c() == n) {flag=1; break;}
%}
if(!flag &&
PDL_IF_BAD($ISGOOD(a()) &&,)
((curind == ns) || $a() '.$op.' cur || PDL_ISNAN_$PPSYM()(cur)))
{cur = $a(); curind = n;}
%}
if (curind != ns) { $c() = curind; }
else { $SETBAD(c()); $PDLSTATESETBAD(c); }
%}',
Doc => <<EOF,
=for refReturns the index of m_size $name elements. As of 2.077, you can specify how many by either passing in an ndarray of the given size (DEPRECATED - will be converted to indx if needed and the input arg will be set to that), or just the size, or a null and the size.
${name}_n_ind(\$pdl, \$out = zeroes(5)); # DEPRECATED
\$out = ${name}_n_ind(\$pdl, 5);
${name}_n_ind(\$pdl, \$out = null, 5);EOF BadDoc => 'Output bad flag is cleared for the output ndarray if sufficient non-bad elements found, else remaining slots in $c() are set bad.
Note that NaNs are considered to be valid values and will "win" over non-NaN; see isfinite and badmask for ways of masking NaNs. ', ); synonym("${name}_n_ind", "${synonym}_n_ind"); } # foreach: $which
pp_def( 'minmaximum', HandleBad => 1, Pars => 'a(n); [o]cmin(); [o] cmax(); indx [o]cmin_ind(); indx [o]cmax_ind();', Code => '$GENERIC() curmin = 0, curmax = 0; /* Handle null ndarray --CED */ PDL_Indx curmin_ind = 0, curmax_ind = 0; int flag = 0; loop(n) %{ PDL_IF_BAD(if ( $ISGOOD(a()) ),) { if ( !flag ) { curmin = curmax = $a(); curmin_ind = curmax_ind = n; flag = 1; } else { if ( $a() < curmin ) { curmin = $a(); curmin_ind = n; } if ( $a() > curmax ) { curmax = $a(); curmax_ind = n; } } } /* ISGOOD */ %} PDL_IF_BAD(if ( !flag ) { /* Handle null ndarray */ $SETBAD(cmin()); $SETBAD(cmin_ind()); $SETBAD(cmax()); $SETBAD(cmax_ind()); $PDLSTATESETBAD(cmin); $PDLSTATESETBAD(cmin_ind); $PDLSTATESETBAD(cmax); $PDLSTATESETBAD(cmax_ind); } else,) { $cmin() = curmin; $cmin_ind() = curmin_ind; $cmax() = curmax; $cmax_ind() = curmax_ind; }', Doc => ' =for ref
Find minimum and maximum and their indices for a given ndarray;
pdl> $x=pdl [[-2,3,4],[1,0,3]]
pdl> ($min, $max, $min_ind, $max_ind)=minmaximum($x)
pdl> p $min, $max, $min_ind, $max_ind
[-2 0] [4 3] [0 1] [2 2]See also "minmax", which clumps the ndarray together.
$name
Return the $text of all elements in an ndarray.
See the documentation for "$func" for more information.
\$x = $name(\$data);This routine handles bad values.
any
Return true if any element in ndarray set
Useful in conditional expressions:
if (any $x>15) { print "some values are greater than 15\n" }See "or" for comments on what happens when all elements in the check are bad.
all
Return true if all elements in ndarray set
Useful in conditional expressions:
if (all $x>15) { print "all values are greater than 15\n" }See "and" for comments on what happens when all elements in the check are bad.
minmax
Returns a list with minimum and maximum values of an ndarray.
($mn, $mx) = minmax($pdl);This routine does not broadcast over the dimensions of $pdl; it returns the minimum and maximum values of the whole ndarray. See "minmaximum" if this is not what is required. The two values are returned as Perl scalars, and therefore ignore whether the values are bad.
pdl> $x = pdl [1,-2,3,5,0]
pdl> ($min, $max) = minmax($x);
pdl> p "$min $max\n";
-2 5$_->[0]pct
Return the specified percentile of all elements in an ndarray. The specified percentile (p) must be between 0.0 and 1.0. When the specified percentile falls between data points, the $_->[1].
\$x = $_->[0]pct(\$data, \$pct);Quicksort a vector into ascending order.
print qsort random(10);Quicksort a vector and return index of elements in ascending order.
$ix = qsorti $x;
print $x->index($ix); # Sorted listSort a list of vectors lexicographically.
The 0th dimension of the source ndarray is dimension in the vector; the 1st dimension is list order. Higher dimensions are broadcasted over.
print qsortvec pdl([[1,2],[0,500],[2,3],[4,2],[3,4],[3,5]]);
[
[ 0 500]
[ 1 2]
[ 2 3]
[ 3 4]
[ 3 5]
[ 4 2]
]
Sort a list of vectors lexicographically, returning the indices of the sorted vectors rather than the sorted list itself.
As with qsortvec, the input PDL should be an NxM array containing M separate N-dimensional vectors. The return value is an integer M-PDL containing the M-indices of original array rows, in sorted order.
As with qsortvec, the zeroth element of the vectors runs slowest in the sorted list.
Additional dimensions are broadcasted over: each plane is sorted separately, so qsortveci may be thought of as a collapse operator of sorts (groan).
AUTHOR
Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu). Contributions by Christian Soeller (c.soeller@auckland.ac.nz) and Karl Glazebrook (kgb@aaoepp.aao.gov.au). All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation under certain conditions. For details, see the file COPYING in the PDL distribution. If this file is separated from the PDL distribution, the copyright notice should be included in the file.