Skip site navigation (1)Skip section navigation (2)

FreeBSD Manual Pages

  
 
  

home | help
Primitive(3)	      User Contributed Perl Documentation	  Primitive(3)

NAME
       PDL::Primitive -	primitive operations for pdl

DESCRIPTION
       This module provides some primitive and useful functions	defined	using
       PDL::PP and able	to use the new indexing	tricks.

       See PDL::Indexing for how to use	indices	creatively.  For explanation
       of the signature	format,	see PDL::PP.

SYNOPSIS
	# Pulls	in PDL::Primitive, among other modules.
	use PDL;

	# Only pull in PDL::Primitive:
	use PDL::Primitive;

FUNCTIONS
   inner
	 Signature: (a(n); b(n); [o]c())

       Inner product over one dimension

	c = sum_i a_i *	b_i

       If "a() * b()" contains only bad	data, "c()" is set bad.	Otherwise
       "c()" will have its bad flag cleared, as	it will	not contain any	bad
       values.

   outer
	 Signature: (a(n); b(m); [o]c(n,m))

       outer product over one dimension

       Naturally, it is	possible to achieve the	effects	of outer product
       simply by threading over	the ""*"" operator but this function is
       provided	for convenience.

       outer processes bad values.  It will set	the bad-value flag of all
       output piddles if the flag is set for any of the	input piddles.

   x
	Signature: (a(i,z), b(x,i),[o]c(x,z))

       Matrix multiplication

       PDL overloads the "x" operator (normally	the repeat operator) for
       matrix multiplication.  The number of columns (size of the 0 dimension)
       in the left-hand	argument must normally equal the number	of rows	(size
       of the 1	dimension) in the right-hand argument.

       Row vectors are represented as (N x 1) two-dimensional PDLs, or you may
       be sloppy and use a one-dimensional PDL.	 Column	vectors	are
       represented as (1 x N) two-dimensional PDLs.

       Threading occurs	in the usual way, but as both the 0 and	1 dimension
       (if present) are	included in the	operation, you must be sure that you
       don't try to thread over	either of those	dims.

       EXAMPLES

       Here are	some simple ways to define vectors and matrices:

	pdl> $r	= pdl(1,2);		   # A row vector
	pdl> $c	= pdl([[3],[4]]);	   # A column vector
	pdl> $c	= pdl(3,4)->(*1);	   # A column vector, using NiceSlice
	pdl> $m	= pdl([[1,2],[3,4]]);	   # A 2x2 matrix

       Now that	we have	a few objects prepared,	here is	how to matrix-multiply
       them:

	pdl> print $r x	$m		   # row x matrix = row
	[
	 [ 7 10]
	]

	pdl> print $m x	$r		   # matrix x row = ERROR
	PDL: Dim mismatch in matmult of	[2x2] x	[2x1]: 2 != 1

	pdl> print $m x	$c		   # matrix x column = column
	[
	 [ 5]
	 [11]
	]

	pdl> print $m x	2		   # Trivial case: scalar mult.
	[
	 [2 4]
	 [6 8]
	]

	pdl> print $r x	$c		   # row x column = scalar
	[
	 [11]
	]

	pdl> print $c x	$r		   # column x row = matrix
	[
	 [3 6]
	 [4 8]
	]

       INTERNALS

       The mechanics of	the multiplication are carried out by the matmult
       method.

   matmult
	 Signature: (a(t,h); b(w,t); [o]c(w,h))

       Matrix multiplication

       Notionally, matrix multiplication $a x $b is equivalent to the
       threading expression

	   $a->dummy(1)->inner($b->xchg(0,1)->dummy(2),$c);

       but for large matrices that breaks CPU cache and	is slow.  Instead,
       matmult calculates its result in	32x32x32 tiles,	to keep	the memory
       footprint within	cache as long as possible on most modern CPUs.

       For usage, see x, a description of the overloaded 'x' operator

       matmult ignores the bad-value flag of the input piddles.	 It will set
       the bad-value flag of all output	piddles	if the flag is set for any of
       the input piddles.

   innerwt
	 Signature: (a(n); b(n); c(n); [o]d())

       Weighted	(i.e. triple) inner product

	d = sum_i a(i) b(i) c(i)

       innerwt processes bad values.  It will set the bad-value	flag of	all
       output piddles if the flag is set for any of the	input piddles.

   inner2
	 Signature: (a(n); b(n,m); c(m); [o]d())

       Inner product of	two vectors and	a matrix

	d = sum_ij a(i)	b(i,j) c(j)

       Note that you should probably not thread	over "a" and "c" since that
       would be	very wasteful. Instead,	you should use a temporary for "b*c".

       inner2 processes	bad values.  It	will set the bad-value flag of all
       output piddles if the flag is set for any of the	input piddles.

   inner2d
	 Signature: (a(n,m); b(n,m); [o]c())

       Inner product over 2 dimensions.

       Equivalent to

	$c = inner($a->clump(2), $b->clump(2))

       inner2d processes bad values.  It will set the bad-value	flag of	all
       output piddles if the flag is set for any of the	input piddles.

   inner2t
	 Signature: (a(j,n); b(n,m); c(m,k); [t]tmp(n,k); [o]d(j,k)))

       Efficient Triple	matrix product "a*b*c"

       Efficiency comes	from by	using the temporary "tmp". This	operation only
       scales as "N**3"	whereas	threading using	inner2 would scale as "N**4".

       The reason for having this routine is that you do not need to have the
       same thread-dimensions for "tmp"	as for the other arguments, which in
       case of large numbers of	matrices makes this much more memory-
       efficient.

       It is hoped that	things like this could be taken	care of	as a kind of
       closures	at some	point.

       inner2t processes bad values.  It will set the bad-value	flag of	all
       output piddles if the flag is set for any of the	input piddles.

   crossp
	 Signature: (a(tri=3); b(tri); [o] c(tri))

       Cross product of	two 3D vectors

       After

	$c = crossp $a,	$b

       the inner product "$c*$a" and "$c*$b" will be zero, i.e.	$c is
       orthogonal to $a	and $b

       crossp does not process bad values.  It will set	the bad-value flag of
       all output piddles if the flag is set for any of	the input piddles.

   norm
	 Signature: (vec(n); [o] norm(n))

       Normalises a vector to unit Euclidean length

       norm processes bad values.  It will set the bad-value flag of all
       output piddles if the flag is set for any of the	input piddles.

   indadd
	 Signature: (a(); indx ind(); [o] sum(m))

       Threaded	Index Add: Add "a" to the "ind"	element	of "sum", i.e:

	sum(ind) += a

       Simple Example:

	 $a = 2;
	 $ind =	3;
	 $sum =	zeroes(10);
	 indadd($a,$ind, $sum);
	 print $sum
	 #Result: ( 2 added to element 3 of $sum)
	 # [0 0	0 2 0 0	0 0 0 0]

       Threaded	Example:

	 $a = pdl( 1,2,3);
	 $ind =	pdl( 1,4,6);
	 $sum =	zeroes(10);
	 indadd($a,$ind, $sum);
	 print $sum."\n";
	 #Result: ( 1, 2, and 3	added to elements 1,4,6	$sum)
	 # [0 1	0 0 2 0	3 0 0 0]

       The routine barfs if any	of the indices are bad.

   conv1d
	 Signature: (a(m); kern(p); [o]b(m); int reflect)

       1D convolution along first dimension

       The m-th	element	of the discrete	convolution of an input	piddle $a of
       size $M,	and a kernel piddle $kern of size $P, is calculated as

				     n = ($P-1)/2
				     ====
				     \
	 ($a conv1d $kern)[m]	=     >	     $a_ext[m -	n] * $kern[n]
				     /
				     ====
				     n = -($P-1)/2

       where $a_ext is either the periodic (or reflected) extension of $a so
       it is equal to $a on " 0..$M-1 "	and equal to the corresponding
       periodic/reflected image	of $a outside that range.

	 $con =	conv1d sequence(10), pdl(-1,0,1);

	 $con =	conv1d sequence(10), pdl(-1,0,1), {Boundary => 'reflect'};

       By default, periodic boundary conditions	are assumed (i.e. wrap
       around).	 Alternatively,	you can	request	reflective boundary conditions
       using the "Boundary" option:

	 {Boundary => 'reflect'} # case	in 'reflect' doesn't matter

       The convolution is performed along the first dimension. To apply	it
       across another dimension	use the	slicing	routines, e.g.

	 $b = $a->mv(2,0)->conv1d($kernel)->mv(0,2); # along third dim

       This function is	useful for threaded filtering of 1D signals.

       Compare also conv2d, convolve, fftconvolve, fftwconv, rfftwconv

       WARNING:	"conv1d" processes bad values in its inputs as the numeric
       value of	"$pdl->badvalue" so it is not recommended for processing pdls
       with bad	values in them unless special care is taken.

       conv1d ignores the bad-value flag of the	input piddles.	It will	set
       the bad-value flag of all output	piddles	if the flag is set for any of
       the input piddles.

   in
	 Signature: (a(); b(n);	[o] c())

       test if a is in the set of values b

	  $goodmsk = $labels->in($goodlabels);
	  print	pdl(3,1,4,6,2)->in(pdl(2,3,3));
	 [1 0 0	0 1]

       "in" is akin to the is an element of of set theory. In principle, PDL
       threading could be used to achieve its functionality by using a
       construct like

	  $msk = ($labels->dummy(0) == $goodlabels)->orover;

       However,	"in" doesn't create a (potentially large) intermediate and is
       generally faster.

       in does not process bad values.	It will	set the	bad-value flag of all
       output piddles if the flag is set for any of the	input piddles.

   uniq
       return all unique elements of a piddle

       The unique elements are returned	in ascending order.

	 PDL> p	pdl(2,2,2,4,0,-1,6,6)->uniq
	 [-1 0 2 4 6]	  # 0 is returned 2nd (sorted order)

	 PDL> p	pdl(2,2,2,4,nan,-1,6,6)->uniq
	 [-1 2 4 6 nan]	  # NaN	value is returned at end

       Note: The returned pdl is 1D; any structure of the input	piddle is
       lost.  "NaN" values are never compare equal to any other	values,	even
       themselves.  As a result, they are always unique. "uniq"	returns	the
       NaN values at the end of	the result piddle.  This follows the Matlab
       usage.

       See uniqind if you need the indices of the unique elements rather than
       the values.

       Bad values are not considered unique by uniq and	are ignored.

	$a=sequence(10);
	$a=$a->setbadif($a%3);
	print $a->uniq;
	[0 3 6 9]

   uniqind
       Return the indices of all unique	elements of a piddle The order is in
       the order of the	values to be consistent	with uniq. "NaN" values	never
       compare equal with any other value and so are always unique.  This
       follows the Matlab usage.

	 PDL> p	pdl(2,2,2,4,0,-1,6,6)->uniqind
	 [5 4 1	3 6]	 # the 0 at index 4 is returned	2nd, but...

	 PDL> p	pdl(2,2,2,4,nan,-1,6,6)->uniqind
	 [5 1 3	6 4]	 # ...the NaN at index 4 is returned at	end

       Note: The returned pdl is 1D; any structure of the input	piddle is
       lost.

       See uniq	if you want the	unique values instead of the indices.

       Bad values are not considered unique by uniqind and are ignored.

   uniqvec
       Return all unique vectors out of	a collection

	 NOTE: If any vectors in the input piddle have NaN values
	 they are returned at the end of the non-NaN ones.  This is
	 because, by definition, NaN values never compare equal	with
	 any other value.

	 NOTE: The current implementation does not sort	the vectors
	 containing NaN	values.

       The unique vectors are returned in lexicographically sorted ascending
       order. The 0th dimension	of the input PDL is treated as a dimensional
       index within each vector, and the 1st and any higher dimensions are
       taken to	run across vectors. The	return value is	always 2D; any
       structure of the	input PDL (beyond using	the 0th	dimension for vector
       index) is lost.

       See also	uniq for a unique list of scalars; and qsortvec	for sorting a
       list of vectors lexicographcally.

       If a vector contains all	bad values, it is ignored as in	uniq.  If some
       of the values are good, it is treated as	a normal vector. For example,
       [1 2 BAD] and [BAD 2 3] could be	returned, but [BAD BAD BAD] could not.
       Vectors containing BAD values will be returned after any	non-NaN	and
       non-BAD containing vectors, followed by the NaN vectors.

   hclip
	 Signature: (a(); b(); [o] c())

       clip (threshold)	$a by $b ($b is	upper bound)

       hclip processes bad values.  It will set	the bad-value flag of all
       output piddles if the flag is set for any of the	input piddles.

   lclip
	 Signature: (a(); b(); [o] c())

       clip (threshold)	$a by $b ($b is	lower bound)

       lclip processes bad values.  It will set	the bad-value flag of all
       output piddles if the flag is set for any of the	input piddles.

   clip
       Clip (threshold)	a piddle by (optional) upper or	lower bounds.

	$b = $a->clip(0,3);
	$c = $a->clip(undef, $x);

       clip handles bad	values since it	is just	a wrapper around hclip and
       lclip.

   clip
	 Signature: (a(); l(); h(); [o]	c())

       info not	available

       clip processes bad values.  It will set the bad-value flag of all
       output piddles if the flag is set for any of the	input piddles.

   wtstat
	 Signature: (a(n); wt(n); avg(); [o]b(); int deg)

       Weighted	statistical moment of given degree

       This calculates a weighted statistic over the vector "a".  The formula
       is

	b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)

       Bad values are ignored in any calculation; $b will only have its	bad
       flag set	if the output contains any bad data.

   statsover
	 Signature: (a(n); w(n); float+	[o]avg(); float+ [o]prms(); int+ [o]median(); int+ [o]min(); int+ [o]max(); float+ [o]adev(); float+ [o]rms())

       Calculate useful	statistics over	a dimension of a piddle

	 ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($piddle, $weights);

       This utility function calculates	various	useful quantities of a piddle.
       These are:

       o  the mean:

	    MEAN = sum (x)/ N

	  with "N" being the number of elements	in x

       o  the population RMS deviation from the	mean:

	    PRMS = sqrt( sum( (x-mean(x))^2 )/(N-1)

	  The population deviation is the best-estimate	of the deviation of
	  the population from which a sample is	drawn.

       o  the median

	  The median is	the 50th percentile data value.	 Median	is found by
	  medover, so WEIGHTING	IS IGNORED FOR THE MEDIAN CALCULATION.

       o  the minimum

       o  the maximum

       o  the average absolute deviation:

	    AADEV = sum( abs(x-mean(x))	)/N

       o  RMS deviation	from the mean:

	    RMS	= sqrt(sum( (x-mean(x))^2 )/N)

	  (also	known as the root-mean-square deviation, or the	square root of
	  the variance)

       This operator is	a projection operator so the calculation will take
       place over the final dimension. Thus if the input is N-dimensional each
       returned	value will be N-1 dimensional, to calculate the	statistics for
       the entire piddle either	use "clump(-1)"	directly on the	piddle or call
       "stats".

       Bad values are simply ignored in	the calculation, effectively reducing
       the sample size.	 If all	data are bad then the output data are marked
       bad.

   stats
       Calculates useful statistics on a piddle

	($mean,$prms,$median,$min,$max,$adev,$rms) = stats($piddle,[$weights]);

       This utility calculates all the most useful quantities in one call.  It
       works the same way as "statsover", except that the quantities are
       calculated considering the entire input PDL as a	single sample, rather
       than as a collection of rows. See "statsover" for definitions of	the
       returned	quantities.

       Bad values are handled; if all input values are bad, then all of	the
       output values are flagged bad.

   histogram
	 Signature: (in(n); int+[o] hist(m); double step; double min; int msize	=> m)

       Calculates a histogram for given	stepsize and minimum.

	$h = histogram($data, $step, $min, $numbins);
	$hist =	zeroes $numbins;  # Put	histogram in existing piddle.
	histogram($data, $hist,	$step, $min, $numbins);

       The histogram will contain $numbins bins	starting from $min, each $step
       wide. The value in each bin is the number of values in $data that lie
       within the bin limits.

       Data below the lower limit is put in the	first bin, and data above the
       upper limit is put in the last bin.

       The output is reset in a	different threadloop so	that you can take a
       histogram of "$a(10,12)"	into "$b(15)" and get the result you want.

       For a higher-level interface, see hist.

	pdl> p histogram(pdl(1,1,2),1,0,3)
	[0 2 1]

       histogram processes bad values.	It will	set the	bad-value flag of all
       output piddles if the flag is set for any of the	input piddles.

   whistogram
	 Signature: (in(n); float+ wt(n);float+[o] hist(m); double step; double	min; int msize => m)

       Calculates a histogram from weighted data for given stepsize and
       minimum.

	$h = whistogram($data, $weights, $step,	$min, $numbins);
	$hist =	zeroes $numbins;  # Put	histogram in existing piddle.
	whistogram($data, $weights, $hist, $step, $min,	$numbins);

       The histogram will contain $numbins bins	starting from $min, each $step
       wide. The value in each bin is the sum of the values in $weights	that
       correspond to values in $data that lie within the bin limits.

       Data below the lower limit is put in the	first bin, and data above the
       upper limit is put in the last bin.

       The output is reset in a	different threadloop so	that you can take a
       histogram of "$a(10,12)"	into "$b(15)" and get the result you want.

	pdl> p whistogram(pdl(1,1,2), pdl(0.1,0.1,0.5),	1, 0, 4)
	[0 0.2 0.5 0]

       whistogram processes bad	values.	 It will set the bad-value flag	of all
       output piddles if the flag is set for any of the	input piddles.

   histogram2d
	 Signature: (ina(n); inb(n); int+[o] hist(ma,mb); double stepa;	double mina; int masize	=> ma;
			    double stepb; double minb; int mbsize => mb;)

       Calculates a 2d histogram.

	$h = histogram2d($datax, $datay, $stepx, $minx,
	      $nbinx, $stepy, $miny, $nbiny);
	$hist =	zeroes $nbinx, $nbiny;	# Put histogram	in existing piddle.
	histogram2d($datax, $datay, $hist, $stepx, $minx,
	      $nbinx, $stepy, $miny, $nbiny);

       The histogram will contain $nbinx x $nbiny bins,	with the lower limits
       of the first one	at "($minx, $miny)", and with bin size "($stepx,
       $stepy)".  The value in each bin	is the number of values	in $datax and
       $datay that lie within the bin limits.

       Data below the lower limit is put in the	first bin, and data above the
       upper limit is put in the last bin.

	pdl> p histogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),1,0,3,1,0,3)
	[
	 [0 0 0]
	 [0 2 2]
	 [0 1 0]
	]

       histogram2d processes bad values.  It will set the bad-value flag of
       all output piddles if the flag is set for any of	the input piddles.

   whistogram2d
	 Signature: (ina(n); inb(n); float+ wt(n);float+[o] hist(ma,mb); double	stepa; double mina; int	masize => ma;
			    double stepb; double minb; int mbsize => mb;)

       Calculates a 2d histogram from weighted data.

	$h = whistogram2d($datax, $datay, $weights,
	      $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
	$hist =	zeroes $nbinx, $nbiny;	# Put histogram	in existing piddle.
	whistogram2d($datax, $datay, $weights, $hist,
	      $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);

       The histogram will contain $nbinx x $nbiny bins,	with the lower limits
       of the first one	at "($minx, $miny)", and with bin size "($stepx,
       $stepy)".  The value in each bin	is the sum of the values in $weights
       that correspond to values in $datax and $datay that lie within the bin
       limits.

       Data below the lower limit is put in the	first bin, and data above the
       upper limit is put in the last bin.

	pdl> p whistogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),pdl(0.1,0.2,0.3,0.4,0.5),1,0,3,1,0,3)
	[
	 [  0	0   0]
	 [  0 0.5 0.9]
	 [  0 0.1   0]
	]

       whistogram2d processes bad values.  It will set the bad-value flag of
       all output piddles if the flag is set for any of	the input piddles.

   fibonacci
	 Signature: ([o]x(n))

       Constructor - a vector with Fibonacci's sequence

       fibonacci does not process bad values.  It will set the bad-value flag
       of all output piddles if	the flag is set	for any	of the input piddles.

   append
	 Signature: (a(n); b(m); [o] c(mn))

       append two or more piddles by concatenating along their first
       dimensions

	$a = ones(2,4,7);
	$b = sequence 5;
	$c = $a->append($b);  #	size of	$c is now (7,4,7) (a jumbo-piddle ;)

       "append"	appends	two piddles along their	first dims. Rest of the
       dimensions must be compatible in	the threading sense. Resulting size of
       first dim is the	sum of the sizes of the	first dims of the two argument
       piddles - ie "n + m".

       Similar functions include glue (below) and cat.

       append does not process bad values.  It will set	the bad-value flag of
       all output piddles if the flag is set for any of	the input piddles.

   glue
	 $c = $a->glue(<dim>,$b,...)

       Glue two	or more	PDLs together along an arbitrary dimension (N-D
       append).

       Sticks $a, $b, and all following	arguments together along the specified
       dimension.  All other dimensions	must be	compatible in the threading
       sense.

       Glue is permissive, in the sense	that every PDL is treated as having an
       infinite	number of trivial dimensions of	order 1	-- so "$a->glue(3,$b)"
       works, even if $a and $b	are only one dimensional.

       If one of the PDLs has no elements, it is ignored.  Likewise, if	one of
       them is actually	the undefined value, it	is treated as if it had	no
       elements.

       If the first parameter is a defined perl	scalar rather than a pdl, then
       it is taken as a	dimension along	which to glue everything else, so you
       can say "$cube =	PDL::glue(3,@image_list);" if you like.

       "glue" is implemented in	pdl, using a combination of xchg and append.
       It should probably be updated (one day) to a pure PP function.

       Similar functions include append	(above)	and cat.

   axisvalues
	 Signature: ([o,nc]a(n))

       Internal	routine

       "axisvalues" is the internal primitive that implements axisvals and
       alters its argument.

       axisvalues does not process bad values.	It will	set the	bad-value flag
       of all output piddles if	the flag is set	for any	of the input piddles.

   random
       Constructor which returns piddle	of random numbers

	$a = random([type], $nx, $ny, $nz,...);
	$a = random $b;

       etc (see	zeroes).

       This is the uniform distribution	between	0 and 1	(assumedly excluding 1
       itself).	The arguments are the same as "zeroes" (q.v.) -	i.e. one can
       specify dimensions, types or give a template.

       You can use the perl function srand to seed the random generator. For
       further details consult Perl's  srand documentation.

   randsym
       Constructor which returns piddle	of random numbers

	$a = randsym([type], $nx, $ny, $nz,...);
	$a = randsym $b;

       etc (see	zeroes).

       This is the uniform distribution	between	0 and 1	(excluding both	0 and
       1, cf random). The arguments are	the same as "zeroes" (q.v.) - i.e. one
       can specify dimensions, types or	give a template.

       You can use the perl function srand to seed the random generator. For
       further details consult Perl's  srand documentation.

   grandom
       Constructor which returns piddle	of Gaussian random numbers

	$a = grandom([type], $nx, $ny, $nz,...);
	$a = grandom $b;

       etc (see	zeroes).

       This is generated using the math	library	routine	"ndtri".

       Mean = 0, Stddev	= 1

       You can use the perl function srand to seed the random generator. For
       further details consult Perl's  srand documentation.

   vsearch
	 Signature: ( vals(); xs(n); [o] indx(); [\%options] )

       Efficiently search for values in	a sorted piddle, returning indices.

	 $idx =	vsearch( $vals,	$x, [\%options]	);
	 vsearch( $vals, $x, $idx, [\%options ]	);

       vsearch performs	a binary search	in the ordered piddle $x, for the
       values from $vals piddle, returning indices into	$x.  What is a
       "match",	and the	meaning	of the returned	indices, are determined	by the
       options.

       The "mode" option indicates which method	of searching to	use, and may
       be one of:

       "sample"
	   invoke vsearch_sample, returning indices appropriate	for sampling
	   within a distribution.

       "insert_leftmost"
	   invoke vsearch_insert_leftmost, returning the left-most possible
	   insertion point which still leaves the piddle sorted.

       "insert_rightmost"
	   invoke vsearch_insert_rightmost, returning the right-most possible
	   insertion point which still leaves the piddle sorted.

       "insert_match"
	   invoke vsearch_match, returning the index of	a matching element,
	   else	-(insertion point + 1)

       "insert_bin_inclusive"
	   invoke vsearch_bin_inclusive, returning an index appropriate	for
	   binning on a	grid where the left bin	edges are inclusive of the
	   bin.	See below for further explanation of the bin.

       "insert_bin_exclusive"
	   invoke vsearch_bin_exclusive, returning an index appropriate	for
	   binning on a	grid where the left bin	edges are exclusive of the
	   bin.	See below for further explanation of the bin.

       The default value of "mode" is "sample".

   vsearch_sample
	 Signature: (vals(); x(n); indx	[o]idx())

       Search for values in a sorted array, return index appropriate for
       sampling	from a distribution

	 $idx =	vsearch_sample($vals, $x);

       $x must be sorted, but may be in	decreasing or increasing order.

       vsearch_sample returns an index I for each value	V of $vals appropriate
       for sampling $vals

       I has the following properties:

       o   if $x is sorted in increasing order

		     V <= x[0]	: I = 0
	     x[0]  < V <= x[-1]	: I s.t. x[I-1]	< V <= x[I]
	     x[-1] < V		: I = $x->nelem	-1

       o   if $x is sorted in decreasing order

		      V	> x[0]	: I = 0
	     x[0]  >= V	> x[-1]	: I s.t. x[I] >= V > x[I+1]
	     x[-1] >= V		: I = $x->nelem	- 1

       If all elements of $x are equal,	I = $x-_nelem -	1.

       If $x contains duplicated elements, I is	the index of the leftmost (by
       position	in array) duplicate if V matches.

       This function is	useful e.g. when you have a list of probabilities for
       events and want to generate indices to events:

	$a = pdl(.01,.86,.93,1); # Barnsley IFS	probabilities cumulatively
	$b = random 20;
	$c = vsearch_sample($b,	$a); # Now, $c will have the appropriate distr.

       It is possible to use the cumusumover function to obtain	cumulative
       probabilities from absolute probabilities.

       needs major (?) work to handles bad values

   vsearch_insert_leftmost
	 Signature: (vals(); x(n); indx	[o]idx())

       Determine the insertion point for values	in a sorted array, inserting
       before duplicates.

	 $idx =	vsearch_insert_leftmost($vals, $x);

       $x must be sorted, but may be in	decreasing or increasing order.

       vsearch_insert_leftmost returns an index	I for each value V of $vals
       equal to	the leftmost position (by index	in array) within $x that V may
       be inserted and still maintain the order	in $x.

       Insertion at index I involves shifting elements I and higher of $x to
       the right by one	and setting the	now empty element at index I to	V.

       I has the following properties:

       o   if $x is sorted in increasing order

		     V <= x[0]	: I = 0
	     x[0]  < V <= x[-1]	: I s.t. x[I-1]	< V <= x[I]
	     x[-1] < V		: I = $x->nelem

       o   if $x is sorted in decreasing order

		      V	>  x[0]	 : I = -1
	     x[0]  >= V	>= x[-1] : I s.t. x[I] >= V > x[I+1]
	     x[-1] >= V		 : I = $x->nelem -1

       If all elements of $x are equal,

	   i = 0

       If $x contains duplicated elements, I is	the index of the leftmost (by
       index in	array) duplicate if V matches.

       needs major (?) work to handles bad values

   vsearch_insert_rightmost
	 Signature: (vals(); x(n); indx	[o]idx())

       Determine the insertion point for values	in a sorted array, inserting
       after duplicates.

	 $idx =	vsearch_insert_rightmost($vals,	$x);

       $x must be sorted, but may be in	decreasing or increasing order.

       vsearch_insert_rightmost	returns	an index I for each value V of $vals
       equal to	the rightmost position (by index in array) within $x that V
       may be inserted and still maintain the order in $x.

       Insertion at index I involves shifting elements I and higher of $x to
       the right by one	and setting the	now empty element at index I to	V.

       I has the following properties:

       o   if $x is sorted in increasing order

		      V	< x[0]	: I = 0
	     x[0]  <= V	< x[-1]	: I s.t. x[I-1]	<= V < x[I]
	     x[-1] <= V		: I = $x->nelem

       o   if $x is sorted in decreasing order

		     V >= x[0]	: I = -1
	     x[0]  > V >= x[-1]	: I s.t. x[I] >= V > x[I+1]
	     x[-1] > V		: I = $x->nelem	-1

       If all elements of $x are equal,

	   i = $x->nelem - 1

       If $x contains duplicated elements, I is	the index of the leftmost (by
       index in	array) duplicate if V matches.

       needs major (?) work to handles bad values

   vsearch_match
	 Signature: (vals(); x(n); indx	[o]idx())

       Match values against a sorted array.

	 $idx =	vsearch_match($vals, $x);

       $x must be sorted, but may be in	decreasing or increasing order.

       vsearch_match returns an	index I	for each value V of $vals.  If V
       matches an element in $x, I is the index	of that	element, otherwise it
       is -( insertion_point + 1 ), where insertion_point is an	index in $x
       where V may be inserted while maintaining the order in $x.  If $x has
       duplicated values, I may	refer to any of	them.

       needs major (?) work to handles bad values

   vsearch_bin_inclusive
	 Signature: (vals(); x(n); indx	[o]idx())

       Determine the index for values in a sorted array	of bins, lower bound
       inclusive.

	 $idx =	vsearch_bin_inclusive($vals, $x);

       $x must be sorted, but may be in	decreasing or increasing order.

       $x represents the edges of contiguous bins, with	the first and last
       elements	representing the outer edges of	the outer bins,	and the	inner
       elements	the shared bin edges.

       The lower bound of a bin	is inclusive to	the bin, its outer bound is
       exclusive to it.	 vsearch_bin_inclusive returns an index	I for each
       value V of $vals

       I has the following properties:

       o   if $x is sorted in increasing order

		      V	< x[0]	: I = -1
	     x[0]  <= V	< x[-1]	: I s.t. x[I] <= V < x[I+1]
	     x[-1] <= V		: I = $x->nelem	- 1

       o   if $x is sorted in decreasing order

		      V	>= x[0]	 : I = 0
	     x[0]  >  V	>= x[-1] : I s.t. x[I+1] > V >=	x[I]
	     x[-1] >  V		 : I = $x->nelem

       If all elements of $x are equal,

	   i = $x->nelem - 1

       If $x contains duplicated elements, I is	the index of the righmost (by
       index in	array) duplicate if V matches.

       needs major (?) work to handles bad values

   vsearch_bin_exclusive
	 Signature: (vals(); x(n); indx	[o]idx())

       Determine the index for values in a sorted array	of bins, lower bound
       exclusive.

	 $idx =	vsearch_bin_exclusive($vals, $x);

       $x must be sorted, but may be in	decreasing or increasing order.

       $x represents the edges of contiguous bins, with	the first and last
       elements	representing the outer edges of	the outer bins,	and the	inner
       elements	the shared bin edges.

       The lower bound of a bin	is exclusive to	the bin, its upper bound is
       inclusive to it.	 vsearch_bin_exclusive returns an index	I for each
       value V of $vals.

       I has the following properties:

       o   if $x is sorted in increasing order

		      V	<= x[0]	 : I = -1
	     x[0]  <  V	<= x[-1] : I s.t. x[I] < V <= x[I+1]
	     x[-1] <  V		 : I = $x->nelem - 1

       o   if $x is sorted in decreasing order

		      V	>  x[0]	 : I = 0
	     x[0]  >= V	>  x[-1] : I s.t. x[I-1] >= V >	x[I]
	     x[-1] >= V		 : I = $x->nelem

       If all elements of $x are equal,

	   i = $x->nelem - 1

       If $x contains duplicated elements, I is	the index of the righmost (by
       index in	array) duplicate if V matches.

       needs major (?) work to handles bad values

   interpolate
	 Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())

       routine for 1D linear interpolation

	( $yi, $err ) =	interpolate($xi, $x, $y)

       Given a set of points "($x,$y)",	use linear interpolation to find the
       values $yi at a set of points $xi.

       "interpolate" uses a binary search to find the suspects,	er...,
       interpolation indices and therefore abscissas (ie $x) have to be
       strictly	ordered	(increasing or decreasing).  For interpolation at lots
       of closely spaced abscissas an approach that uses the last index	found
       as a start for the next search can be faster (compare Numerical Recipes
       "hunt" routine).	Feel free to implement that on top of the binary
       search if you like. For out of bounds values it just does a linear
       extrapolation and sets the corresponding	element	of $err	to 1, which is
       otherwise 0.

       See also	interpol, which	uses the same routine, differing only in the
       handling	of extrapolation - an error message is printed rather than
       returning an error piddle.

       needs major (?) work to handles bad values

   interpol
	Signature: (xi(); x(n);	y(n); [o] yi())

       routine for 1D linear interpolation

	$yi = interpol($xi, $x,	$y)

       "interpol" uses the same	search method as interpolate, hence $x must be
       strictly	ordered	(either	increasing or decreasing).  The	difference
       occurs in the handling of out-of-bounds values; here an error message
       is printed.

   interpND
       Interpolate values from an N-D piddle, with switchable method

	 $source = 10*xvals(10,10) + yvals(10,10);
	 $index	= pdl([[2.2,3.5],[4.1,5.0]],[[6.0,7.4],[8,9]]);
	 print $source->interpND( $index );

       InterpND	acts like indexND, collapsing $index by	lookup into $source;
       but it does interpolation rather	than direct sampling.  The
       interpolation method and	boundary condition are switchable via an
       options hash.

       By default, linear or sample interpolation is used, with	constant value
       outside the boundaries of the source pdl.  No dataflow occurs, because
       in general the output is	computed rather	than indexed.

       All the interpolation methods treat the pixels as value-centered, so
       the "sample" method will	return "$a->(0)" for coordinate	values on the
       set [-0.5,0.5), and all methods will return "$a->(1)" for a coordinate
       value of	exactly	1.

       Recognized options:

       method
	  Values can be:

	  o  0,	s, sample, Sample (default for integer source types)

	     The nearest value is taken. Pixels	are regarded as	centered on
	     their respective integer coordinates (no offset from the linear
	     case).

	  o  1,	l, linear, Linear (default for floating	point source types)

	     The values	are N-linearly interpolated from an N-dimensional cube
	     of	size 2.

	  o  3,	c, cube, cubic,	Cubic

	     The values	are interpolated using a local cubic fit to the	data.
	     The fit is	constrained to match the original data and its
	     derivative	at the data points.  The second	derivative of the fit
	     is	not continuous at the data points.  Multidimensional datasets
	     are interpolated by the successive-collapse method.

	     (Note that	the constraint on the first derivative causes a	small
	     amount of ringing around sudden features such as step functions).

	  o  f,	fft, fourier, Fourier

	     The source	is Fourier transformed,	and the	interpolated values
	     are explicitly calculated from the	coefficients.  The boundary
	     condition option is ignored -- periodic boundaries	are imposed.

	     If	you pass in the	option "fft", and it is	a list (ARRAY) ref,
	     then it is	a stash	for the	magnitude and phase of the source FFT.
	     If	the list has two elements then they are	taken as already
	     computed; otherwise they are calculated and put in	the stash.

       b, bound, boundary, Boundary
	  This option is passed	unmodified into	indexND, which is used as the
	  indexing engine for the interpolation.  Some current allowed values
	  are 'extend',	'periodic', 'truncate',	and 'mirror' (default is
	  'truncate').

       bad
	  contains the fill value used for 'truncate' boundary.	 (default 0)

       fft
	  An array ref whose associated	list is	used to	stash the FFT of the
	  source data, for the FFT method.

   one2nd
       Converts	a one dimensional index	piddle to a set	of ND coordinates

	@coords=one2nd($a, $indices)

       returns an array	of piddles containing the ND indexes corresponding to
       the one dimensional list	indices. The indices are assumed to correspond
       to array	$a clumped using "clump(-1)". This routine is used in the old
       vector form of whichND, but is useful on	its own	occasionally.

       Returned	piddles	have the indx datatype.	 $indices can have values
       larger than "$a->nelem" but negative values in $indices will not	give
       the answer you expect.

	pdl> $a=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$a->clump(-1)
	pdl> $maxind=maximum_ind($c); p	$maxind;
	6
	pdl> print one2nd($a, maximum_ind($c))
	0 1 1
	pdl> p $a->at(0,1,1)
	3

   which
	 Signature: (mask(n); indx [o] inds(m))

       Returns indices of non-zero values from a 1-D PDL

	$i = which($mask);

       returns a pdl with indices for all those	elements that are nonzero in
       the mask. Note that the returned	indices	will be	1D. If you feed	in a
       multidimensional	mask, it will be flattened before the indices are
       calculated.  See	also whichND for multidimensional masks.

       If you want to index into the original mask or a	similar	piddle with
       output from "which", remember to	flatten	it before calling index:

	 $data = random	5, 5;
	 $idx =	which $data > 0.5; # $idx is now 1D
	 $bigsum = $data->flat->index($idx)->sum;  # flatten before indexing

       Compare also where for similar functionality.

       SEE ALSO:

       which_both returns separately the indices of both zero and nonzero
       values in the mask.

       where returns associated	values from a data PDL,	rather than indices
       into the	mask PDL.

       whichND returns N-D indices into	a multidimensional PDL.

	pdl> $x	= sequence(10);	p $x
	[0 1 2 3 4 5 6 7 8 9]
	pdl> $indx = which($x>6); p $indx
	[7 8 9]

       which processes bad values.  It will set	the bad-value flag of all
       output piddles if the flag is set for any of the	input piddles.

   which_both
	 Signature: (mask(n); indx [o] inds(m);	indx [o]notinds(q))

       Returns indices of zero and nonzero values in a mask PDL

	($i, $c_i) = which_both($mask);

       This works just as which, but the complement of $i will be in $c_i.

	pdl> $x	= sequence(10);	p $x
	[0 1 2 3 4 5 6 7 8 9]
	pdl> ($small, $big) = which_both ($x >=	5); p "$small\n	$big"
	[5 6 7 8 9]
	[0 1 2 3 4]

       which_both processes bad	values.	 It will set the bad-value flag	of all
       output piddles if the flag is set for any of the	input piddles.

   where
       Use a mask to select values from	one or more data PDLs

       "where" accepts one or more data	piddles	and a mask piddle.  It returns
       a list of output	piddles, corresponding to the input data piddles.
       Each output piddle is a 1-dimensional list of values in its
       corresponding data piddle. The values are drawn from locations where
       the mask	is nonzero.

       The output PDLs are still connected to the original data	PDLs, for the
       purpose of dataflow.

       "where" combines	the functionality of which and index into a single
       operation.

       BUGS:

       While "where" works OK for most N-dimensional cases, it does not	thread
       properly	over (for example) the (N+1)th dimension in data that is
       compared	to an N-dimensional mask.  Use "whereND" for that.

	$i = $x->where($x+5 > 0); # $i contains	those elements of $x
				  # where mask ($x+5 > 0) is 1
	$i .= -5;  # Set those elements	(of $x)	to -5. Together, these
		   # commands clamp $x to a maximum of -5.

       It is also possible to use the same mask	for several piddles with the
       same call:

	($i,$j,$k) = where($x,$y,$z, $x+5>0);

       Note: $i	is always 1-D, even if $x is >1-D.

       WARNING:	The first argument (the	values)	and the	second argument	(the
       mask) currently have to have the	exact same dimensions (or horrible
       things happen). You *cannot* thread over	a smaller mask,	for example.

   whereND
       "where" with support for	ND masks and threading

       "whereND" accepts one or	more data piddles and a	mask piddle.  It
       returns a list of output	piddles, corresponding to the input data
       piddles.	 The values are	drawn from locations where the mask is
       nonzero.

       "whereND" differs from "where" in that the mask dimensionality is
       preserved which allows for proper threading of the selection operation
       over higher dimensions.

       As with "where" the output PDLs are still connected to the original
       data PDLs, for the purpose of dataflow.

	 $sdata	= whereND $data, $mask
	 ($s1, $s2, ..., $sn) =	whereND	$d1, $d2, ..., $dn, $mask

	 where

	   $data is M dimensional
	   $mask is N <	M dimensional
	   dims($data) 1..N == dims($mask) 1..N
	   with	threading over N+1 to M	dimensions

	 $data	 = sequence(4,3,2);   #	example	data array
	 $mask4	 = (random(4)>0.5);   #	example	1-D mask array,	has $n4	true values
	 $mask43 = (random(4,3)>0.5); #	example	2-D mask array,	has $n43 true values
	 $sdat4	 = whereND $data, $mask4;   # $sdat4 is	a [$n4,3,2] pdl
	 $sdat43 = whereND $data, $mask43;  # $sdat43 is a [$n43,2] pdl

       Just as with "where", you can use the returned value in an assignment.
       That means that both of these examples are valid:

	 # Used	to create a new	slice stored in	$sdat4:
	 $sdat4	= $data->whereND($mask4);
	 $sdat4	.= 0;
	 # Used	in lvalue context:
	 $data->whereND($mask4)	.= 0;

   whichND
       Return the coordinates of non-zero values in a mask.

       WhichND returns the N-dimensional coordinates of	each nonzero value in
       a mask PDL with any number of dimensions.  The returned values arrive
       as an array-of-vectors suitable for use in indexND or range.

	$coords	= whichND($mask);

       returns a PDL containing	the coordinates	of the elements	that are non-
       zero in $mask, suitable for use in indexND.  The	0th dimension contains
       the full	coordinate listing of each point; the 1st dimension lists all
       the points.  For	example, if $mask has rank 4 and 100 matching
       elements, then $coords has dimension 4x100.

       If no such elements exist, then whichND returns a structured empty PDL:
       an Nx0 PDL that contains	no values (but matches,	threading-wise,	with
       the vectors that	would be produced if such elements existed).

       DEPRECATED BEHAVIOR IN LIST CONTEXT:

       whichND once delivered different	values in list context than in scalar
       context,	for historical reasons.	 In list context, it returned the
       coordinates transposed, as a collection of 1-PDLs (one per dimension)
       in a list.  This	usage is deprecated in PDL 2.4.10, and will cause a
       warning to be issued every time it is encountered.  To avoid the
       warning,	you can	set the	global variable	"$PDL::whichND"	to 's' to get
       scalar behavior in all contexts,	or to 'l' to get list behavior in list
       context.

       In later	versions of PDL, the deprecated	behavior will disappear.
       Deprecated list context whichND expressions can be replaced with:

	   @list = $a->whichND->mv(0,-1)->dog;

       SEE ALSO:

       which finds coordinates of nonzero values in a 1-D mask.

       where extracts values from a data PDL that are associated with nonzero
       values in a mask	PDL.

	pdl> $a=sequence(10,10,3,4)
	pdl> ($x, $y, $z, $w)=whichND($a == 203); p $x,	$y, $z,	$w
	[3] [0]	[2] [0]
	pdl> print $a->at(list(cat($x,$y,$z,$w)))
	203

   setops
       Implements simple set operations	like union and intersection

	  Usage: $set =	setops($a, <OPERATOR>, $b);

       The operator can	be "OR", "XOR" or "AND". This is then applied to $a
       viewed as a set and $b viewed as	a set. Set theory says that a set may
       not have	two or more identical elements,	but setops takes care of this
       for you,	so "$a=pdl(1,1,2)" is OK. The functioning is as	follows:

       "OR"
	   The resulting vector	will contain the elements that are either in
	   $a or in $b or both.	This is	the union in set operation terms

       "XOR"
	   The resulting vector	will contain the elements that are either in
	   $a or $b, but not in	both. This is

		Union($a, $b) -	Intersection($a, $b)

	   in set operation terms.

       "AND"
	   The resulting vector	will contain the intersection of $a and	$b, so
	   the elements	that are in both $a and	$b. Note that for convenience
	   this	operation is also aliased to intersect.

       It should be emphasized that these routines are used when one or	both
       of the sets $a, $b are hard to calculate	or that	you get	from a
       separate	subroutine.

       Finally IDL users might be familiar with	Craig Markwardt's
       "cmset_op.pro" routine which has	inspired this routine although it was
       written independently However the present routine has a few less
       options (but see	the examples)

       You will	very often use these functions on an index vector, so that is
       what we will show here. We will in fact something slightly silly. First
       we will find all	squares	that are also cubes below 10000.

       Create a	sequence vector:

	 pdl> $x = sequence(10000)

       Find all	odd and	even elements:

	 pdl> ($even, $odd) = which_both( ($x %	2) == 0)

       Find all	squares

	 pdl> $squares=	which(ceil(sqrt($x)) ==	floor(sqrt($x)))

       Find all	cubes (being careful with roundoff error!)

	 pdl> $cubes= which(ceil($x**(1.0/3.0))	== floor($x**(1.0/3.0)+1e-6))

       Then find all squares that are cubes:

	 pdl> $both = setops($squares, 'AND', $cubes)

       And print these (assumes	that "PDL::NiceSlice" is loaded!)

	 pdl> p	$x($both)
	  [0 1 64 729 4096]

       Then find all numbers that are either cubes or squares, but not both:

	 pdl> $cube_xor_square = setops($squares, 'XOR', $cubes)

	 pdl> p	$cube_xor_square->nelem()
	  112

       So there	are a total of 112 of these!

       Finally find all	odd squares:

	 pdl> $odd_squares = setops($squares, 'AND', $odd)

       Another common occurrence is to want to get all objects that are	in $a
       and in the complement of	$b. But	it is almost always best to create the
       complement explicitly since the universe	that both are taken from is
       not known. Thus use which_both if possible to keep track	of
       complements.

       If this is impossible the best approach is to make a temporary:

       This creates an index vector the	size of	the universe of	the sets and
       set all elements	in $b to 0

	 pdl> $tmp = ones($n_universe);	$tmp($b) .= 0;

       This then finds the complement of $b

	 pdl> $C_b = which($tmp	== 1);

       and this	does the final selection:

	 pdl> $set = setops($a,	'AND', $C_b)

   intersect
       Calculate the intersection of two piddles

	  Usage: $set =	intersect($a, $b);

       This routine is merely a	simple interface to setops. See	that for more
       information

       Find all	numbers	less that 100 that are of the form 2*y and 3*x

	pdl> $x=sequence(100)
	pdl> $factor2 =	which( ($x % 2)	== 0)
	pdl> $factor3 =	which( ($x % 3)	== 0)
	pdl> $ii=intersect($factor2, $factor3)
	pdl> p $x($ii)
	[0 6 12	18 24 30 36 42 48 54 60	66 72 78 84 90 96]

AUTHOR
       Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu).
       Contributions by	Christian Soeller (c.soeller@auckland.ac.nz), Karl
       Glazebrook (kgb@aaoepp.aao.gov.au), Craig DeForest
       (deforest@boulder.swri.edu) and Jarle Brinchmann	(jarle@astro.up.pt)
       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.

       Updated for CPAN	viewing	compatibility by David Mertens.

perl v5.24.1			  2017-07-02			  Primitive(3)

NAME | DESCRIPTION | SYNOPSIS | FUNCTIONS | AUTHOR

Want to link to this manual page? Use this URL:
<https://www.freebsd.org/cgi/man.cgi?query=PDL::Primitive&sektion=3&manpath=FreeBSD+12.1-RELEASE+and+Ports>

home | help