Skip to content

Commit

Permalink
test 08 in idiomatic PDL
Browse files Browse the repository at this point in the history
  • Loading branch information
mohawk2 committed Jan 15, 2023
1 parent 638379a commit f481834
Showing 1 changed file with 25 additions and 39 deletions.
64 changes: 25 additions & 39 deletions t/x08.pl
Original file line number Diff line number Diff line change
Expand Up @@ -98,54 +98,40 @@ sub cmap1_init {
my $x = (sequence (XPTS) - int(XPTS / 2)) / int(XPTS / 2);
$x *= 1.5
if $rosen;
my $xx = $x->dummy(1,YPTS);

my $y = (sequence (YPTS) - int(YPTS / 2)) / int(YPTS / 2);
$y += 0.5
if $rosen;
my $yy = $y->dummy(0,XPTS);

my $z = zeroes (XPTS, YPTS);
my ($i, $j);
for ($i = 0; $i < XPTS; $i++) {
my $xx = $x->index ($i);
for ($j = 0; $j < YPTS; $j++) {
my $yy = $y->index ($j);
my $zz;
if ($rosen) {
$zz = (1 - $xx) ** 2 + 100 * ($yy - ($xx ** 2)) ** 2;
# The log argument may be zero for just the right grid.
if ($zz > 0.) {
$zz = log ($zz);
} else {
$zz = -5.; # -MAXFLOAT would mess-up up the scale
}
}
else {
my $r = sqrt ($xx * $xx + $yy * $yy);
$zz = exp (-$r * $r) * cos (2.0 * pi * $r);
}
$z->index ($i)->index ($j) .= $zz;
}
my $z;
if ($rosen) {
$z = (1 - $xx) ** 2 + 100 * ($yy - ($xx ** 2)) ** 2;
# The log argument may be zero for just the right grid.
$z = log ($z);
}

sub mymin { $_[0] < $_[1] ? $_[0] : $_[1] }
sub mymax { $_[0] > $_[1] ? $_[0] : $_[1] }

else {
my $r = sqrt ($xx * $xx + $yy * $yy);
$z = exp (-$r * $r) * cos (2.0 * pi * $r);
}
$z->inplace->setnonfinitetobad;
$z->inplace->setbadtoval(-5); # -MAXFLOAT would mess-up up the scale

my (@indexymin, @indexymax);
my $i = sequence( XPTS );
my $square_root = sqrt( 1. - hclip( ( ( $i - $x0 ) / $a ) ** 2, 1 ) );
# Add 0.5 to find nearest integer and therefore preserve symmetry
# with regard to lower and upper bound of y range.
my $indexymin = lclip( 0.5 + $y0 - $b * $square_root, 0 )->indx;
# indexymax calculated with the convention that it is 1
# greater than highest valid index.
my $indexymax = hclip( 1 + ( 0.5 + $y0 + $b * $square_root ), YPTS )->indx;
my $zlimited = zeroes (XPTS, YPTS);
my (@indexymin, @indexymax, @zlimited);
for my $i ( $indexxmin..$indexxmax-1 ) {
my $square_root = sqrt( 1. - mymin( 1, ( ( $i - $x0 ) / $a ) ** 2 ) );
# Add 0.5 to find nearest integer and therefore preserve symmetry
# with regard to lower and upper bound of y range.
$indexymin[$i] = mymax( 0, ( 0.5 + $y0 - $b * $square_root ) );
# indexymax calculated with the convention that it is 1
# greater than highest valid index.
$indexymax[$i] = mymin( YPTS, 1 + ( 0.5 + $y0 + $b * $square_root ) );
for my $j ($indexymin[$i]..$indexymax[$i]-1) {
$zlimited->index($i)->index($j) .= $z->index($i)->index($j);
}
my $j = [ $indexymin->at($i), $indexymax->at($i) ];
$zlimited->index($i)->slice($j) .= $z->index($i)->slice($j);
}
my $indexymin = pdl \@indexymin;
my $indexymax = pdl \@indexymax;

my $zmin = min ($z);
my $zmax = max ($z);
Expand Down

0 comments on commit f481834

Please sign in to comment.