3
# check for bad value support
5
my $bvalflag = $PDL::Config{WITH_BADVAL} || 0;
10
#error "You must have a working RAND_MAX! Something's wrong with your include files"
15
pp_addpm({At=>'Top'},<<'EOD');
22
PDL::Primitive - primitive operations for pdl
26
This module provides some primitive and useful functions defined
27
using PDL::PP and able to use the new indexing tricks.
29
See L<PDL::Indexing|PDL::Indexing> for how to use indices creatively.
30
For explanation of the signature format, see L<PDL::PP|PDL::PP>.
40
pp_addpm({At=>'Bot'},<<'EOD');
44
Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu). Contributions
45
by Christian Soeller (c.soeller@auckland.ac.nz) and Karl Glazebrook
46
(kgb@aaoepp.aao.gov.au).
47
All rights reserved. There is no warranty. You are allowed
48
to redistribute this software / documentation under certain
49
conditions. For details, see the file COPYING in the PDL
50
distribution. If this file is separated from the PDL distribution,
51
the copyright notice should be included in the file.
58
################################################################
59
# a whole bunch of quite basic functions for inner, outer
60
# and matrix products (operations that are not normally
61
# available via operator overloading)
62
################################################################
67
Pars => 'a(n); b(n); [o]c();',
70
loop(n) %{ tmp += $a() * $b(); %}
76
if ( $ISGOOD(a()) && $ISGOOD(b()) ) { tmp += $a() * $b(); flag = 1; }
78
if ( flag ) { $c() = tmp; }
79
else { $SETBAD(c()); $PDLSTATESETBAD(c); }',
80
CopyBadStatusCode => '',
84
Inner product over one dimension
90
'If C<a() * b()> contains only bad data,
91
C<c()> is set bad. Otherwise C<c()> will have its bad flag cleared,
92
as it will not contain any bad values.',
99
Pars => 'a(n); b(m); [o]c(n,m);',
106
if ( $ISBAD(a()) || $ISBAD(b()) ) {
115
outer product over one dimension
117
Naturally, it is possible to achieve the effects of outer
118
product simply by threading over the "C<*>"
119
operator but this function is provided for convenience.
121
'); # pp_def( outer )
129
Signature: (a(x,y),b(y,z),[o]c(x,z))
133
Matrix multiplication
135
We peruse the inner product to define matrix multiplication
136
via a threaded inner product
141
barf "Invalid number of arguments for matmult" if $#_ < 1;
143
while ($a->getndims < 2) {$a = $a->dummy(-1)} # promote if necessary
144
while ($b->getndims < 2) {$b = $b->dummy(-1)}
145
if(!defined $c) {$c = PDL->nullcreate($a)}
146
$a->dummy(1)->inner($b->xchg(0,1)->dummy(2),$c);
150
*matmult = \&PDL::matmult;
154
pp_add_exported('', 'matmult');
159
Pars => 'a(n); b(n); c(n); [o]d();',
163
tmp += $a() * $b() * $c();
171
if ( $ISGOOD(a()) && $ISGOOD(b()) && $ISGOOD(c()) ) {
172
tmp += $a() * $b() * $c();
176
if ( flag ) { $d() = tmp; }
177
else { $SETBAD(d()); }',
181
Weighted (i.e. triple) inner product
183
d = sum_i a(i) b(i) c(i)
191
Pars => 'a(n); b(n,m); c(m); [o]d();',
195
tmp += $a() * $b() * $c();
202
if ( $ISGOOD(a()) && $ISGOOD(b()) && $ISGOOD(c()) ) {
203
tmp += $a() * $b() * $c();
207
if ( flag ) { $d() = tmp; }
208
else { $SETBAD(d()); }',
212
Inner product of two vectors and a matrix
214
d = sum_ij a(i) b(i,j) c(j)
216
Note that you should probably not thread over C<a> and C<c> since that would be
217
very wasteful. Instead, you should use a temporary for C<b*c>.
226
Pars => 'a(n,m); b(n,m); [o]c();',
237
if ( $ISGOOD(a()) && $ISGOOD(b()) ) {
242
if ( flag ) { $c() = tmp; }
243
else { $SETBAD(c()); }',
247
Inner product over 2 dimensions.
251
$c = inner($a->clump(2), $b->clump(2))
259
Pars => 'a(j,n); b(n,m); c(m,k); [t]tmp(n,k); [o]d(j,k));',
271
tmp1 += $a() * $tmp();
280
if ( $ISGOOD(b()) && $ISGOOD(c()) ) {
285
if ( flag ) { $tmp() = tmp0; }
286
else { $SETBAD(tmp()); }
292
if ( $ISGOOD(a()) && $ISGOOD(tmp()) ) {
293
tmp1 += $a() * $tmp();
297
if ( flag ) { $d() = tmp1; }
298
else { $SETBAD(d()); }
303
Efficient Triple matrix product C<a*b*c>
305
Efficiency comes from by using the temporary C<tmp>. This operation only
306
scales as C<N**3> whereas threading using L<inner2|/inner2> would scale
309
The reason for having this routine is that you do not need to
310
have the same thread-dimensions for C<tmp> as for the other arguments,
311
which in case of large numbers of matrices makes this much more
314
It is hoped that things like this could be taken care of as a kind of
315
closures at some point.
318
); # pp_def inner2t()
321
# a helper function for the cross product definition
323
"\$c(tri => $_[0]) = \$a(tri => $_[1])*\$b(tri => $_[2]) -
324
\$a(tri => $_[2])*\$b(tri => $_[1]);"
331
Cross product of two 3D vectors
339
the inner product C<$c*$a> and C<$c*$b> will be zero, i.e. C<$c> is
340
orthogonal to C<$a> and C<$b>
344
Pars => 'a(tri=3); b(tri); [o] c(tri)',
353
Pars => 'vec(n); [o] norm(n)',
354
Doc => 'Normalises a vector to unit Euclidean length',
357
loop(n) %{ sum += $vec()*$vec(); %}
360
loop(n) %{ $norm() = $vec()/sum; %}
362
loop(n) %{ $norm() = $vec(); %}
368
sum += $vec()*$vec();
375
if ( $ISBAD(vec()) ) { $SETBAD(norm()); }
376
else { $norm() = $vec()/sum; }
380
if ( $ISBAD(vec()) ) { $SETBAD(norm()); }
381
else { $norm() = $vec(); }
391
# this one was motivated by the need to compute
392
# the circular mean efficiently
393
# without it could not be done efficiently or without
394
# creating large intermediates (check pdl-porters for
396
# see PDL::ImageND for info about the circ_mean function
400
Pars => 'a(); int ind(); [o] sum(m)',
402
'register int foo = $ind();
403
if( foo<0 || foo>=$SIZE(m) ) {
404
barf("PDL::indadd: invalid index");
406
$sum(m => foo) += $a();',
408
'register int foo = $ind();
409
if( $ISBADVAR(foo,ind) || foo<0 || foo>=$SIZE(m) ) {
410
barf("PDL::indadd: invalid index");
412
if ( $ISBAD(a()) ) { $SETBAD(sum(m => foo)); }
413
else { $sum(m => foo) += $a(); }',
415
'The routine barfs if any of the indices are bad.',
419
Threaded Index Add: Add C<a> to the C<ind> element of C<sum>, i.e:
430
indadd($a,$ind, $sum);
432
#Result: ( 2 added to element 3 of $sum)
433
# [0 0 0 2 0 0 0 0 0 0]
440
indadd($a,$ind, $sum);
442
#Result: ( 1, 2, and 3 added to elements 1,4,6 $sum)
443
# [0 1 0 0 2 0 3 0 0 0]
449
# useful for threaded 1D filters
451
/* Fast Modulus with proper negative behaviour */
452
#define REALMOD(a,b) while ((a)>=(b)) (a) -= (b); while ((a)<0) (a) += (b);
458
1d convolution along first dimension
462
$con = conv1d sequence(10), pdl(-1,0,1), {Boundary => 'reflect'};
464
By default, periodic boundary conditions are assumed (i.e. wrap around).
465
Alternatively, you can request reflective boundary conditions using
466
the C<Boundary> option:
468
{Boundary => 'reflect'} # case in 'reflect' doesn't matter
470
The convolution is performed along the first dimension. To apply it across
471
another dimension use the slicing routines, e.g.
473
$b = $a->mv(2,0)->conv1d($kernel)->mv(0,2); # along third dim
475
This function is useful for threaded filtering of 1D signals.
477
Compare also L<conv2d|PDL::Image2D/conv2d>, L<convolve|PDL::ImageND/convolve>,
478
L<fftconvolve|PDL::FFT/fftconvolve>, L<fftwconv|PDL::FFTW/fftwconv>,
479
L<rfftwconv|PDL::FFTW/rfftwconv>
482
Pars => 'a(m); kern(p); [o]b(m);',
483
OtherPars => 'int reflect;',
487
my $opt = pop @_ if ref($_[$#_]) eq \'HASH\';
488
die \'Usage: conv1d( a(m), kern(p), [o]b(m), {Options} )\'
491
my $c = $#_ == 2 ? $_[2] : PDL->null;
492
&PDL::_conv1d_int($a,$kern,$c,
493
!(defined $opt && exists $$opt{Boundary}) ? 0 :
494
lc $$opt{Boundary} eq "reflect");
502
int reflect = $COMP(reflect);
503
int m_size = $COMP(__m_size);
504
int p_size = $COMP(__p_size);
507
for(i=0; i<m_size; i++) {
509
for(i1=0; i1<p_size; i1++) {
513
if (reflect && i2>=m_size)
514
i2 = m_size-(i2-m_size+1);
517
tmp += $a(m=>i2) * $kern(p=>i1);
524
# this can be achieved by
525
# ($a->dummy(0) == $b)->orover
526
# but this one avoids a larger intermediate and potentially shortcuts
528
Pars => 'a(); b(n); [o] c()',
530
loop(n) %{ if ($a() == $b()) {$c() = 1; break;} %}',
534
test if a is in the set of values b
538
$goodmsk = $labels->in($goodlabels);
539
print pdl(4,3,1)->in(pdl(2,3,3));
542
C<in> is akin to the I<is an element of> of set theory. In priciple,
543
PDL threading could be used to achieve its functionality by using a
546
$msk = ($labels->dummy(0) == $goodlabels)->orover;
548
However, C<in> doesn't create a (potentially large) intermediate
549
and is generally faster.
554
pp_add_exported '', 'uniq';
561
return all unique elements of a piddle
563
The unique elements are returned in ascending order.
567
print pdl(2,2,2,4,0,-1,6,6)->uniq;
570
Note: The returned pdl is 1D; any structure of the input
576
# return unique elements of array
577
# find as jumps in the sorted array
578
# flattens in the process
580
use PDL::Core 'barf';
582
barf "need at least one element" if $arr->nelem == 0;
583
my $srt = $arr->clump(-1)->qsort;
584
my $uniq = $srt->where($srt != $srt->rotate(-1));
585
# make sure we return something if there is only one value
586
return $uniq->nelem == 0 ? $arr->pdl($arr->type,[$arr->clump(-1)->at(0)]) :
592
#####################################################################
594
#####################################################################
602
my $name = $opt->[0];
607
Pars => 'a(); b(); [o] c()',
609
'$c() = ($a() '.$op.' $b()) ? $b() : $a();',
611
'if ( $ISBAD(a()) || $ISBAD(b()) ) {
614
$c() = ($a() '.$op.' $b()) ? $b() : $a();
616
Doc => 'clip C<$a> by C<$b> (C<$b> is '.
617
($name eq 'hclip' ? 'upper' : 'lower').' bound)',
622
if (\$a->is_inplace) {
623
\$a->set_inplace(0); \$c = \$a;
624
} elsif (\$#_ > 1) {\$c=\$_[2]} else {\$c=PDL->nullcreate(\$a)}
625
&PDL::_${name}_int(\$a,\$b,\$c);
633
pp_add_exported('', 'clip');
641
Clip a piddle by (optional) upper or lower bounds.
646
$c = $a->clip(undef, $x);
654
clip handles bad values since it is just a
655
wrapper around L<hclip|/hclip> and
667
my $d; if($a->is_inplace) {$a->set_inplace(0); $d = $a}
668
elsif($#_ > 2) {$d=$_[3]} else {$d = PDL->nullcreate($a)}
670
&PDL::_lclip_int($a,$b,$d);
672
&PDL::_hclip_int($d,$c,$d);
674
} elsif(defined $c) {
675
&PDL::_hclip_int($a,$c,$d);
682
############################################################
683
# elementary statistics and histograms
684
############################################################
689
Pars => 'a(n); wt(n); avg(); [o]b();',
690
OtherPars => 'int deg',
699
for(i=0; i<$COMP(deg); i++)
701
statsum += $wt() * (tmp - $avg());
703
$b() = statsum / wtsum;',
709
if ( $ISGOOD(wt()) && $ISGOOD(a()) && $ISGOOD(avg()) ) {
714
for(i=0; i<$COMP(deg); i++)
716
statsum += $wt() * (tmp - $avg());
720
if ( flag ) { $b() = statsum / wtsum; }
721
else { $SETBAD(b()); $PDLSTATESETBAD(b); }',
722
CopyBadStatusCode => '',
726
Weighted statistical moment of given degree
728
This calculates a weighted statistic over the vector C<a>.
731
b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)
735
'Bad values are ignored in any calculation; C<$b> will only
736
have its bad flag set if the output contains any bad data.',
743
Pars => 'a(n); w(n); int+ [o]avg(); int+ [o]rms(); int+ [o]min(); int+ [o]max(); int+ [o]adev()',
745
'$GENERIC(avg) tmp = 0;
746
$GENERIC(rms) tmp1 = 0;
747
$GENERIC(adev) tmp2 = 0;
748
$GENERIC(rms) diff = 0;
749
$GENERIC(min) curmin, curmax;
750
$GENERIC(w) norm = 0;
754
if (!n) { curmin = $a(); curmax = $a();}
757
} else if ($a() > curmax) {
764
/* Calculate the RMS using the corrected two-pass algorithm */
767
diff = ($a()*$w() - $avg());
772
$rms() = sqrt ( (tmp - tmp1*tmp1/norm)/((norm - ($GENERIC(rms)) 1)) );
773
$adev() = sqrt( tmp2/ norm );',
775
'$GENERIC(avg) tmp = 0;
776
$GENERIC(rms) tmp1 = 0;
777
$GENERIC(adev) tmp2 = 0;
778
$GENERIC(rms) diff = 0;
779
$GENERIC(min) curmin, curmax;
780
$GENERIC(w) norm = 0;
783
/* perhaps should check w() for bad values too ? */
784
if ( $ISGOOD(a()) ) {
787
if (!flag) { curmin = $a(); curmax = $a(); flag=1; }
790
} else if ($a() > curmax) {
795
/* have at least one valid point if flag == 1 */
803
diff = ($a()*$w() - $avg());
809
$rms() = sqrt ( (tmp - tmp1*tmp1/norm)/( (norm - ($GENERIC(rms)) 1)));
810
$adev() = sqrt ( tmp2 / norm);
821
barf(\'Usage: ($mean,[$rms, $median, $min, $max, $adev]) = statsover($data,[$weights])\') if $#_>1;
822
my ($data, $weights) = @_;
823
$weights = $data->ones() if !defined($weights);
825
my $median = $data->medover();
826
my $mean = PDL->nullcreate($data);
827
my $rms = PDL->nullcreate($data);
828
my $min = PDL->nullcreate($data);
829
my $max = PDL->nullcreate($data);
830
my $adev = PDL->nullcreate($data);
831
&PDL::_statsover_int($data, $weights, $mean, $rms, $min, $max, $adev);
833
return $mean unless wantarray;
834
return ($mean, $rms, $median, $min, $max, $adev);
842
Calculate useful statistics over a dimension of a piddle
846
($mean, $rms, $median, $min, $max, $adev) = statover($piddle, $weights);
848
This utility function calculates various useful
849
quantities of a piddle. These are the mean:
853
with C<N> being the number of elements in x, the root mean
854
square deviation from the mean, RMS, given as,
856
RMS = sqrt(sum( (x-mean(x))^2 )/(N-1));
858
Note the use of C<N-1> which for almost all cases should be
859
the right normalisation factor. The routine also returns
860
the median, minimum and maximum of the piddle as well as
861
the mean absolute deviation, defined as:
863
ADEV = sqrt(sum( abs(x-mean(x)) )/N)
865
note here that we use the mean and not the median. This could
866
possibly be changed in future versions of the code.
868
This operator is a projection operator so the calculation
869
will take place over the final dimension. Thus if the input
870
is N-dimensional each returned value will be N-1 dimensional,
871
to calculate the statistics for the entire piddle either
872
use C<clump(-1)> directly on the piddle or call C<stats>.
878
Bad values are simply ignored in the calculation and if all data are
879
bad the results will also be bad.
883
pp_add_exported('','stats');
890
Calculates useful statistics on a piddle
894
($mean,$rms,$median,$min,$max) = stats($piddle,[$weights]);
896
This utility calculates all the most useful
897
quantities in one call.
899
B<Note:> The RMS value that this function returns in the RMS
900
deviation from the mean, also known as the population standard-
909
The return values if all elements are bad is currently poorly
917
*stats = \&PDL::stats;
919
barf('Usage: ($mean,[$rms]) = stats($data,[$weights])') if $#_>1;
920
my ($data,$weights) = @_;
926
pp_addpm(<<'!NO!SUBS!');
927
my $npts = $data->ngood;
930
pp_addpm(<<'!NO!SUBS!');
931
my $npts = $data->nelem;
938
$mean = 0.0; $rms = 0.0;
940
$mean = $data->dsum / $npts;
941
$rms = sqrt( ((($data-$mean)**2 )->dsum) / $npts );
945
my $wtsum = $weights->dsum;
946
if ( $wtsum == 0.0 ) {
947
$mean = 0.0; $rms = 0.0;
949
$mean = (($weights*$data)->dsum) / $wtsum;
950
$rms = sqrt( ( ( $weights*(($data-$mean)**2) )->dsum ) / $wtsum );
953
my ($median,$min,$max) = ($data->median,$data->min,$data->max);
954
print "Mean = $mean, RMS = $rms, Median = $median\n".
955
"Min = $min, Max = $max\n" if $PDL::verbose;
956
return $mean unless wantarray;
957
return ($mean,$rms,$median,$min,$max);
962
{Name => 'histogram',
968
Doc3 => "number of\n",
969
Doc4 => "\nUse L<hist|PDL::Basic/hist> instead for a high-level interface.\n",
970
Doc5 => "histogram(pdl(1,1,2),1,0,3)\n [0 2 1]"
972
{Name => 'whistogram',
973
WeightPar => 'float+ wt(n);',
974
HistType => 'float+',
975
HistOp => '+= $wt()',
976
Doc1 => " from weighted data",
977
Doc2 => "\$weights, ",
978
Doc3 => "sum of the values in C<\$weights>\nthat correspond to ",
980
Doc5 => "whistogram(pdl(1,1,2), pdl(0.1,0.1,0.5), 1, 0, 4)\n [0 0.2 0.5 0]"
985
Pars => 'in(n); '.$_->{WeightPar}.$_->{HistType}. '[o] hist(m)',
987
OtherPars => 'double step; double min; int msize => m',
991
register int maxj = $SIZE(m)-1;
992
register double min = $COMP(min);
993
register double step = $COMP(step);
995
loop(m) %{ $hist() = 0; %}
999
j = (int) (($in()-min)/step);
1001
if (j > maxj) j = maxj;
1002
($hist(m => j))'.$_->{HistOp}.';
1007
register int maxj = $SIZE(m)-1;
1008
register double min = $COMP(min);
1009
register double step = $COMP(step);
1011
loop(m) %{ $hist() = 0; %}
1015
if ( $ISGOOD(in()) ) {
1016
j = (int) (($in()-min)/step);
1018
if (j > maxj) j = maxj;
1019
($hist(m => j))'.$_->{HistOp}.';
1027
Calculates a histogram$_->{Doc1} for given stepsize and minimum.
1031
\$h = $_->{Name}(\$data, $_->{Doc2}\$step, \$min, \$numbins);
1032
\$hist = zeroes \$numbins; # Put histogram in existing piddle.
1033
$_->{Name}(\$data, $_->{Doc2}\$hist, \$step, \$min, \$numbins);
1035
The histogram will contain C<\$numbins> bins starting from C<\$min>, each
1036
C<\$step> wide. The value in each bin is the $_->{Doc3}values in C<\$data> that lie within the bin limits.
1038
Data below the lower limit is put in the first bin, and data above the
1039
upper limit is put in the last bin.
1041
The output is reset in a different threadloop so that you
1042
can take a histogram of C<\$a(10,12)> into C<\$b(15)> and get the result
1048
perldl> p $_->{Doc5}
1054
{Name => 'histogram2d',
1060
Doc3 => "number of\n",
1061
Doc5 => "histogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),1,0,3,1,0,3)
1069
{Name => 'whistogram2d',
1070
WeightPar => 'float+ wt(n);',
1071
HistType => 'float+',
1072
HistOp => '+= $wt()',
1073
Doc1 => " from weighted data",
1074
Doc2 => " \$weights,",
1075
Doc3 => "sum of the values in\nC<\$weights> that correspond to ",
1076
Doc5 => "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)
1087
Pars => 'ina(n); inb(n); '.$_->{WeightPar}.$_->{HistType}. '[o] hist(ma,mb)',
1088
# set outdim by Par!
1089
OtherPars => 'double stepa; double mina; int masize => ma;
1090
double stepb; double minb; int mbsize => mb;',
1093
'register int ja,jb;
1094
register int maxja = $SIZE(ma)-1;
1095
register int maxjb = $SIZE(mb)-1;
1096
register double mina = $COMP(mina);
1097
register double minb = $COMP(minb);
1098
register double stepa = $COMP(stepa);
1099
register double stepb = $COMP(stepb);
1101
loop(ma,mb) %{ $hist() = 0; %}
1105
ja = (int) (($ina()-mina)/stepa);
1106
jb = (int) (($inb()-minb)/stepb);
1108
if (ja > maxja) ja = maxja;
1110
if (jb > maxjb) jb = maxjb;
1111
($hist(ma => ja,mb => jb))'.$_->{HistOp}.';
1116
'register int ja,jb;
1117
register int maxja = $SIZE(ma)-1;
1118
register int maxjb = $SIZE(mb)-1;
1119
register double mina = $COMP(mina);
1120
register double minb = $COMP(minb);
1121
register double stepa = $COMP(stepa);
1122
register double stepb = $COMP(stepb);
1124
loop(ma,mb) %{ $hist() = 0; %}
1128
if ( $ISGOOD(ina()) && $ISGOOD(inb()) ) {
1129
ja = (int) (($ina()-mina)/stepa);
1130
jb = (int) (($inb()-minb)/stepb);
1132
if (ja > maxja) ja = maxja;
1134
if (jb > maxjb) jb = maxjb;
1135
($hist(ma => ja,mb => jb))'.$_->{HistOp}.';
1144
Calculates a 2d histogram$_->{Doc1}.
1148
\$h = $_->{Name}(\$datax, \$datay,$_->{Doc2}
1149
\$stepx, \$minx, \$nbinx, \$stepy, \$miny, \$nbiny);
1150
\$hist = zeroes \$nbinx, \$nbiny; # Put histogram in existing piddle.
1151
$_->{Name}(\$datax, \$datay,$_->{Doc2} \$hist,
1152
\$stepx, \$minx, \$nbinx, \$stepy, \$miny, \$nbiny);
1154
The histogram will contain C<\$nbinx> x C<\$nbiny> bins, with the lower
1155
limits of the first one at C<(\$minx, \$miny)>, and with bin size
1156
C<(\$stepx, \$stepy)>.
1157
The value in each bin is the $_->{Doc3}values in C<\$datax> and C<\$datay> that lie within the bin limits.
1159
Data below the lower limit is put in the first bin, and data above the
1160
upper limit is put in the last bin.
1164
perldl> p $_->{Doc5}
1170
###########################################################
1171
# a number of constructors: fibonacci, append, axisvalues &
1173
###########################################################
1177
Doc=>'Constructor - a vector with Fibonacci\'s sequence',
1180
sub fibonacci { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->fibonacci : PDL->fibonacci(@_) }
1183
my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inplace;
1184
&PDL::_fibonacci_int($x->clump(-1));
1204
Pars => 'a(n); b(m); [o] c(mn)',
1205
# note that ideally we want to say '$SIZE(mn) = $SIZE(m)+$SIZE(n);'
1206
# but that requires placing RedoDimsParsedCode *after* assignment of
1207
# childdims to $SIZE(XXX)!!! XXXXXmake that workXXXXX
1209
pdl * dpdla = $PDL(a);
1210
pdl * dpdlb = $PDL(b);
1211
$SIZE(mn) = (dpdla->ndims > 0 ? dpdla->dims[0] : 1) +
1212
(dpdlb->ndims > 0 ? dpdlb->dims[0] : 1);
1214
Code => 'register PDL_Long mnp;
1215
PDL_Long ns = $SIZE(n);
1217
loop(n) %{ $c(mn => n) = $a(); %}
1218
loop(m) %{ mnp = m+ns; $c(mn => mnp) = $b(); %}
1223
append two piddles by concantening along their respective first dimensions
1229
$c = $a->append($b); # size of $c is now (7,4,7) (a jumbo-piddle ;)
1231
C<append> appends two piddles along their first dims. Rest of the dimensions
1232
must be compatible in the threading sense. Resulting size of first dim is
1233
the sum of the sizes of the first dims of the two argument piddles -
1239
pp_def( 'axisvalues',
1240
Pars => '[o,nc]a(n)',
1241
Code => 'loop(n) %{ $a() = n; %}',
1247
C<axisvalues> is the internal primitive that implements
1248
L<axisvals|PDL::Basic/axisvals>
1249
and alters its argument.
1252
); # pp_def: axisvalues
1261
Constructor which returns piddle of random numbers
1265
$a = random([type], $nx, $ny, $nz,...);
1268
etc (see L<zeroes|PDL::Core/zeroes>).
1270
This is the uniform distribution between 0 and 1 (assumedly
1271
excluding 1 itself). The arguments are the same as C<zeroes>
1272
(q.v.) - i.e. one can specify dimensions, types or give
1275
You can use the perl function L<srand|perlfunc/srand> to seed the random
1276
generator. For further details consult Perl's L<srand|perlfunc/srand>
1283
Constructor which returns piddle of random numbers
1287
$a = randsym([type], $nx, $ny, $nz,...);
1290
etc (see L<zeroes|PDL::Core/zeroes>).
1292
This is the uniform distribution between 0 and 1 (excluding both 0 and
1293
1, cf L<random|/random>). The arguments are the same as C<zeroes> (q.v.) -
1294
i.e. one can specify dimensions, types or give a template.
1296
You can use the perl function L<srand|perlfunc/srand> to seed the random
1297
generator. For further details consult Perl's L<srand|perlfunc/srand>
1306
#define Drand01() (((double)rand()) / (RAND_MAX+1.0))
1316
'$a() = Drand01();',
1319
sub random { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->random : PDL->random(@_) }
1322
my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inplace;
1323
&PDL::_random_int($x);
1335
do tmp = Drand01(); while (tmp == 0.0); /* 0 < tmp < 1 */
1339
sub randsym { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->randsym : PDL->randsym(@_) }
1342
my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inplace;
1343
&PDL::_randsym_int($x);
1355
Constructor which returns piddle of Gaussian random numbers
1359
$a = grandom([type], $nx, $ny, $nz,...);
1362
etc (see L<zeroes|PDL::Core/zeroes>).
1364
This is generated using the math library routine C<ndtri>.
1366
Mean = 0, Stddev = 1
1369
You can use the perl function L<srand|perlfunc/srand> to seed the random
1370
generator. For further details consult Perl's L<srand|perlfunc/srand>
1375
sub grandom { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->grandom : PDL->grandom(@_) }
1378
my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inplace;
1379
use PDL::Math 'ndtri';
1380
$x .= ndtri(randsym($x));
1386
pp_add_exported('','grandom');
1388
###############################################################
1389
# routines somehow related to interpolation
1390
###############################################################
1392
# The last x is ignored...
1395
BadDoc => 'needs major (?) work to handles bad values',
1396
Pars => 'i(); x(n); int [o]ip()',
1397
GenericTypes => ['F','D'], # too restrictive ?
1398
Code => 'int carp=0;
1400
long n1 = $SIZE(n)-1;
1401
long jl=-1, jh=n1, m;
1402
int up = ($x(n => n1) > $x(n => 0));
1404
while (jh-jl > 1) /* binary search */
1407
if ($i() > $x(n => m) == up)
1414
} else if (jl == n1) {
1415
if ($i() != $x(n => n1)) carp = 1;
1422
if (carp) warn("some values had to be extrapolated");
1427
routine for searching 1D values i.e. step-function interpolation.
1431
$inds = vsearch($vals, $xs);
1433
Returns for each value of C<$vals> the index of the least larger member
1434
of C<$xs> (which need to be in increasing order). If the value is larger
1435
than any member of C<$xs>, the index to the last element of C<$xs> is
1440
This function is useful e.g. when you have a list of probabilities
1441
for events and want to generate indices to events:
1443
$a = pdl(.01,.86,.93,1); # Barnsley IFS probabilities cumulatively
1445
$c = vsearch($b, $a); # Now, $c will have the appropriate distr.
1447
It is possible to use the L<cumusumover|/cumusumover> function to obtain
1448
cumulative probabilities from absolute probabilities.
1452
pp_def('interpolate',
1454
BadDoc => 'needs major (?) work to handles bad values',
1455
Pars => 'xi(); x(n); y(n); [o] yi(); int [o] err()',
1456
GenericTypes => ['F','D'], # too restrictive ?
1461
int up = ($x(n => n1) > $x(n => 0));
1469
while (jh-jl > 1) /* binary search */
1472
if ($xi() > $x(n => m) == up)
1478
if ($xi() != $x(n => 0)) carp = 1;
1480
} else if (jh == n) {
1481
if ($xi() != $x(n => n1)) carp = 1;
1485
if ((d = $x(n => jh)-$x(n => jl)) == 0)
1486
barf("identical abscissas");
1487
d = ($x(n => jh)-$xi())/d;
1488
$yi() = d*$y(n => jl) + (1-d)*$y(n => jh);
1495
routine for 1D linear interpolation
1499
( $yi, $err ) = interpolate($xi, $x, $y)
1501
Given a set of points C<($x,$y)>, use linear interpolation
1502
to find the values C<$yi> at a set of points C<$xi>.
1504
C<interpolate> uses a binary search to find the suspects, er...,
1505
interpolation indices and therefore abscissas (ie C<$x>)
1506
have to be I<strictly> ordered (increasing or decreasing).
1507
For interpolation at lots of
1508
closely spaced abscissas an approach that uses the last index found as
1509
a start for the next search can be faster (compare Numerical Recipes
1510
C<hunt> routine). Feel free to implement that on top of the binary
1511
search if you like. For out of bounds values it just does a linear
1512
extrapolation and sets the corresponding element of C<$err> to 1,
1513
which is otherwise 0.
1515
See also L<interpol|/interpol>, which uses the same routine,
1516
differing only in the handling of extrapolation - an error message
1517
is printed rather than returning an error piddle.
1522
pp_add_exported('', 'interpol');
1529
Signature: (xi(); x(n); y(n); [o] yi())
1533
routine for 1D linear interpolation
1537
$yi = interpol($xi, $x, $y)
1539
C<interpol> uses the same search method as L<interpolate|/interpolate>,
1540
hence C<$x> must be I<strictly> ordered (either increasing or decreasing).
1541
The difference occurs in the handling of out-of-bounds values; here
1542
an error message is printed.
1546
# kept in for backwards compatability
1547
sub interpol ($$$;$) {
1553
if ( $#_ == 0 ) { $yi = $_[0]; }
1554
else { $yi = PDL->null; }
1556
interpolate( $xi, $x, $y, $yi, my $err = PDL->null );
1557
print "some values had to be extrapolated\n"
1560
return $yi if $#_ == -1;
1563
*PDL::interpol = \&interpol;
1568
##############################################################
1569
# things related to indexing: one2nd, which, where
1570
##############################################################
1572
pp_add_exported("", 'one2nd');
1579
Converts a one dimensional index piddle to a set of ND coordinates
1583
@coords=one2nd($a, $indices)
1585
returns an array of piddles containing the ND indexes corresponding to
1586
the one dimensional list indices. The indices are assumed to correspond
1587
to array C<$a> clumped using C<clump(-1)>. This routine is used in
1588
L<whichND|/whichND>,
1589
but is useful on its own occasionally.
1593
perldl> $a=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$a->clump(-1)
1594
perldl> $maxind=maximum_ind($c); p $maxind;
1596
perldl> print one2nd($a, maximum_ind($c))
1598
perldl> p $a->at(0,1,1)
1603
*one2nd = \&PDL::one2nd;
1605
barf "Usage: one2nd $array $indices\n" if $#_ != 1;
1607
my @dimension=$a->dims;
1610
foreach (@dimension) {
1611
$index[$count++]=$ind % $_;
1619
my $doc_which = <<'EOD';
1623
Returns piddle of indices of non-zero values.
1629
returns a pdl with indices for all those elements that are
1630
nonzero in the mask. Note that the returned indices will be 1D. If you want
1631
to index into the original mask or a similar piddle remember to flatten
1632
it before calling index:
1634
$data = random 5, 5;
1635
$idx = which $data > 0.5; # $idx is now 1D
1636
$bigsum = $data->flat->index($idx)->sum; # flatten before indexing
1638
Compare also L<where|/where> for similar functionality.
1640
If you want to return both the indices of non-zero values and the
1641
complement, use the function L<which_both|/which_both>.
1645
perldl> $x = sequence(10); p $x
1646
[0 1 2 3 4 5 6 7 8 9]
1647
perldl> $indx = which($x>6); p $indx
1652
my $doc_which_both = <<'EOD';
1656
Returns piddle of indices of non-zero values and their complement
1660
($i, $c_i) = which_both($mask);
1662
This works just as L<which|/which>, but the complement of C<$i> will be in
1667
perldl> $x = sequence(10); p $x
1668
[0 1 2 3 4 5 6 7 8 9]
1669
perldl> ($small, $big) = which_both ($x >= 5); p "$small\n $big"
1677
Pars => 'mask(n); int [o] inds(m);',
1678
Variables => 'int dm=0;',
1680
Autosize => '$SIZE(m) = sum;',
1683
sub which { my ($this,$out) = @_;
1684
$this = $this->flat;
1685
$out = $this->nullcreate unless defined $out;
1686
PDL::_which_int($this,$out);
1689
*PDL::which = \&which;
1692
{Name => 'which_both',
1693
Pars => 'mask(n); int [o] inds(m); int [o]notinds(q)',
1694
Variables => 'int dm=0; int dm2=0;',
1695
Elseclause => "else { \n \$notinds(q => dm2)=n; \n dm2++;\n }",
1696
Autosize => '$SIZE(m) = sum;'."\n".' $SIZE(q) = dpdl->dims[0]-sum;',
1697
Doc => $doc_which_both,
1699
sub which_both { my ($this,$outi,$outni) = @_;
1700
$this = $this->flat;
1701
$outi = $this->nullcreate unless defined $outi;
1702
$outni = $this->nullcreate unless defined $outni;
1703
PDL::_which_both_int($this,$outi,$outni);
1704
return wantarray ? ($outi,$outni) : $outi;
1706
*PDL::which_both = \&which_both;
1714
PMCode => $_->{PMCode},
1715
Code => $_->{Variables} .
1720
}'.$_->{Elseclause} . "\n".
1722
# the next one is currently a dirty hack
1723
# this will probably break once dataflow is enabled again
1724
# *unless* we have made sure that mask is physical by now!!!
1727
/* not sure if this is necessary */
1728
pdl * dpdl = $PDL(mask);
1729
$GENERIC() *m_datap = (($GENERIC() *)(PDL_REPRP(dpdl)));
1730
PDL_Long inc = PDL_REPRINC(dpdl,0);
1731
PDL_Long offs = PDL_REPROFFS(dpdl);
1734
if (dpdl->ndims != 1)
1735
barf("dimflag currently works only with 1D pdls");
1737
for (i=0; i<dpdl->dims[0]; i++) {
1738
if ( *(m_datap+inc*i+offs)) sum++;
1741
'.$_->{Autosize} . '
1742
/* printf("RedoDimsCode: setting dim m to %ld\n",sum); */'
1751
Returns indices to non-zero values or those values from another piddle.
1755
$i = $x->where($x+5 > 0); # $i contains elements of $x
1756
# where mask ($x+5 > 0) is 1
1758
Note: C<$i> is always 1-D, even if C<$x> is >1-D. The first argument
1759
(the values) and the second argument (the mask) currently have to have
1760
the same initial dimensions (or horrible things happen).
1762
It is also possible to use the same mask for several piddles with
1765
($i,$j,$k) = where($x,$y,$z, $x+5>0);
1770
barf "Usage: where( $pdl1, ..., $pdlN, $mask )\n" if $#_ == 0;
1773
my($data,$mask) = @_;
1774
$data = $_[0]->clump(-1) if $_[0]->getndims>1;
1775
$mask = $_[1]->clump(-1) if $_[0]->getndims>1;
1776
return $data->index($mask->which());
1778
if($_[-1]->getndims > 1) {
1779
my $mask = $_[-1]->clump(-1)->which;
1780
return map {$_->clump(-1)->index($mask)} @_[0..$#_-1];
1782
my $mask = $_[-1]->which;
1783
return map {$_->index($mask)} @_[0..$#_-1];
1787
*where = \&PDL::where;
1792
pp_add_exported("", 'where');
1799
Returns the coordinates for non-zero values
1803
@coords=whichND($mask);
1805
returns an array of piddles containing the coordinates of the elements
1806
that are non-zero in C<$mask>.
1810
perldl> $a=sequence(10,10,3,4)
1811
perldl> ($x, $y, $z, $w)=whichND($a == 203); p $x, $y, $z, $w
1813
perldl> print $a->at(list(cat($x,$y,$z,$w)))
1818
*whichND = \&PDL::whichND;
1821
my $ind=($mask->clump(-1))->which;
1823
return $mask->one2nd($ind);
1829
pp_add_exported("", 'whichND');
added
removed
Basic/Primitive/primitive.pd
