Skip to content

Commit

Permalink
Slatec, Proj4 tests replace tapprox with is_pdl - #34
Browse files Browse the repository at this point in the history
  • Loading branch information
@mohawk2
mohawk2 committed Oct 30, 2024
1 parent 3a664ee commit 21979bb
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 70 deletions.
93 changes: 36 additions & 57 deletions Libtmp/Slatec/t/slatec.t
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,45 +2,36 @@ use strict;
use warnings;
use PDL::LiteF;
use Test::More;
use Test::PDL;
use PDL::Slatec;
use PDL::MatrixOps qw(identity);

kill 'INT', $$ if $ENV{UNDER_DEBUGGER}; # Useful for debugging.

sub tapprox {
my($x,$y,$c,$d) = @_;
$c = abs($x-$y);
$d = max($c);
# print "APR: $x,$y,$c,$d;\n";
$d < 0.001;
}

my $mat = pdl [1,0.1],[0.1,2];

my ($eigvals,$eigvecs) = eigsys($mat);

## print STDERR $eigvecs,$eigvals,"\n";

ok(tapprox($eigvals,pdl(0.9901,2.009)));
ok(!tapprox($eigvals,pdl(0.99,2.5)));

ok(tapprox($eigvecs,pdl([0.995,-0.0985],[0.0985,0.995])));
is_pdl $eigvals,float(0.9901,2.009), {atol=>1e-3};
is_pdl $eigvecs,float([0.995,-0.0985],[0.0985,0.995]), {atol=>1e-3};

$mat = pdl [2,3],[4,5];

my $inv = matinv($mat);

my $uni=scalar $mat x $inv;
ok(tapprox($uni,identity(2)));
is_pdl $uni,identity(2);

eval {matinv(identity(2)->dummy(-1,2))};
is $@, '', 'matinv can broadcast';

my $det = $mat->det;
my $deti = $inv->det;

ok(tapprox($det,-2));
ok(tapprox($deti,-0.5));
is_pdl $det, pdl(-2);
is_pdl $deti, pdl(-0.5);

# Now do the polynomial fitting tests

Expand Down Expand Up @@ -133,38 +124,34 @@ $xx = pdl(12,4,6.25,1.5); # Ask for multiple positions at once

($yfit, $yp) = polyvalue($ndeg, $nder, $xx, $a1);

## print STDERR "At $xx is $yfit and $yp\n";

# Simple test of expected value
ok(int($yfit->at(1)) == 15);

# test the PCHIP stuff

## Test: chim/chic
#
$x = float( 3 .. 10 );
my $f = $x*$x*$x + 425.42352;
my $answer = 3.0*$x*$x;
my $answer = 3*$x*$x;

my ( $d, $err ) = chim( float($x), float($f) );

ok(($err->getndims==0) & ($err->sum == 0));
is_pdl $err, longlong 0;

# don't check the first and last elements, as expect the
# error to be largest there
# value of 5% comes from tests on linux and solaris machines
ok(all( slice( abs(($d - $answer)/$answer), '1:-2' ) < 0.05 ) );
is_pdl +(map $_->slice('1:-2'), $d, $answer), {rtol=>0.05};

# compare the results of chic
my $d2 = $f->zeroes;
chic( pdl([0, 0]), pdl([0, 0]), 1, $x, $f, $d2, my $err2=null );
ok(($err2->getndims==0) & ($err2->sum == 0));
ok(all( abs($d2 - $d) < 0.02 ) );
is_pdl $err2, longlong 0;
is_pdl $d2, $d->double, {atol=>2e-2};

## Test: chsp
chsp( pdl([0, 0]), pdl([0, 0]), $x, $f, my $d3=null, my $err3=null );
ok(($err3->getndims==0) & ($err3->sum == 0));
ok(all( abs($d3 - $d) < 2 ) );
is_pdl $err3, longlong 0;
is_pdl $d3, $d->double, {atol=>2};

## Test: chfd/chfe
#
Expand All @@ -175,33 +162,27 @@ my ( $fe, $de );
ok(($err->getndims==0) & ($err->sum == 0));

$answer = $xe*$xe*$xe + 425.42352;
ok(all( abs(($fe - $answer)/$answer) < 1.0e-5 ) );
is_pdl $fe, $answer, {rtol=>1e-5};

$answer = 3.0*$xe*$xe;
ok(all( abs(($de - $answer)/$answer) < 0.02 ) );
is_pdl $de, $answer, {rtol=>2e-2};

( $fe, $err ) = chfe( $x, $f, $d, 1, $xe );

ok(($err->getndims==0) & ($err->sum == 0));

$answer = $xe*$xe*$xe + 425.42352;
ok(all( abs(($fe - $answer)/$answer) < 1.0e-5 ) );
is_pdl $fe, $xe*$xe*$xe + 425.42352, {rtol=>1e-3};
is_pdl $err, longlong 0;

# Test: chcm
#
$x = float( 1, 2, 3, 5, 6, 7 );
$f = float( 1, 2, 3, 4, 3, 4 );
my $ans = long( 1, 1, 1, -1, 1, 2 );
my $ans = longlong( 1, 1, 1, -1, 1, 2 );

( $d, $err ) = chim($x, $f);
ok(($err->getndims==0) & ($err->sum == 2)); # 2 switches in monotonicity

my $ismon;
( $ismon, $err ) = chcm($x, $f, $d, 1);
is_pdl $err, longlong 2;

ok(($err->getndims==0) & ($err->sum == 0));
ok($ismon->get_datatype == longlong->enum) or diag $ismon->type;
ok(tapprox($ismon,$ans)) or diag $ismon, "\nexpected:\n", $ans;
( my $ismon, $err ) = chcm($x, $f, $d, 1);
is_pdl $err, longlong 0;
is_pdl $ismon,$ans;

## Test: chia
#
Expand All @@ -211,17 +192,15 @@ $f = $x * $x;

$ans = pdl( 9.0**3, (8.0**3-1.0**3) ) / 3.0;
( my $int, $err ) = chia($x, $f, $d, my $skip=zeroes(2), pdl(0.0,1.0), pdl(9.0,8.0));
ok(all($err == 0));
ok(all( abs($int-$ans) < 0.04 ) );
is_pdl $err, longlong 0,0;
is_pdl $int, $ans, {atol=>4e-2};

my $hi = pdl( $x->at(9), $x->at(7) );
my $lo = pdl( $x->at(0), $x->at(1) );
$ans = ($hi**3 - $lo**3) / 3;
( $int, $err ) = chid( $x, $f, $d, $skip=zeroes(2), pdl(0,1), pdl(9,7) );
ok(all($err == 0));
ok(all( abs($int-$ans) < 0.06 ) );
## print STDERR "int=$int; ans=$ans; int-ans=".($int-$ans)."\n";
## print STDERR "ref ans=".(ref $ans)."\n";
is_pdl $err, longlong 0,0;
is_pdl $int, $ans, {atol=>6e-2};

## Test: chbs - note, only tests that it runs successfully
my $nknots = 0;
Expand All @@ -231,7 +210,7 @@ my $ndim = PDL->null;
my $kord = PDL->null;
$err = PDL->null;
chbs( $x, $f, $d, 0, $nknots, $t, $bcoef, $ndim, $kord, $err );
ok(all($err == 0));
is_pdl $err, longlong 0;

## Test: bvalu - note, only tests that it runs successfully
my $x_slice = $x->slice('0:-2'); # because calling with last value is out of range
Expand All @@ -244,21 +223,21 @@ gefa(my $lu=$A->copy, my $ipiv=null, my $info=null);
gesl($lu, $ipiv, $x=$B->transpose->copy, 1); # 1 = do transpose because Fortran
$x = $x->inplace->transpose;
my $got = $A x $x;
ok tapprox $got, $B or diag "got: $got";
is_pdl $got, $B;

{
my $pa = pdl(float,1,-1,1,-1); # even number
my ($az, $x, $y) = PDL::Slatec::fft($pa);
ok all approx $az, 0;
ok all approx $x, pdl "[0 1 0 0]";
ok all approx $y, pdl "[0 0 0 0]";
ok all approx PDL::Slatec::rfft($az, $x, $y), $pa;
is_pdl $az, float 0;
is_pdl $x, float "[0 1 0 0]";
is_pdl $y, float "[0 0 0 0]";
is_pdl PDL::Slatec::rfft($az, $x, $y), $pa;
$pa = pdl(float,1,-1,1,-1,1); # odd number
($az, $x, $y) = PDL::Slatec::fft($pa);
ok all approx $az, 0.2;
ok all approx $x, pdl "[0.4 0.4 0 0 0]";
ok all approx $y, pdl "[-0.2906170 -1.231073 0 0 0]";
ok all approx PDL::Slatec::rfft($az, $x, $y), $pa;
is_pdl $az, float 0.2;
is_pdl $x, float "[0.4 0.4 0 0 0]";
is_pdl $y, float "[-0.2906170 -1.231073 0 0 0]";
is_pdl PDL::Slatec::rfft($az, $x, $y), $pa;
}

done_testing;
17 changes: 4 additions & 13 deletions Libtmp/Transform-Proj4/t/gis_proj.t
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ use strict;
use warnings;
use PDL::LiteF;
use Test::More;
use Test::PDL;
BEGIN {
diag "ENV $_ = ", explain $ENV{$_}
for qw(LD_LIBRARY_PATH DYLD_LIBRARY_PATH LDFLAGS CFLAGS CXXFLAGS LD_RUN_PATH
Expand All @@ -11,16 +12,6 @@ use PDL::Transform::Proj4;
use Alien::proj;
diag "Alien::proj version $Alien::proj::VERSION";

sub tapprox
{
my $x = shift;
my $y = shift;
my $d = abs($x - $y);
my $res = all($d < 1.0e-5);
diag "got:$x\nexpected:$y" if !$res;
$res;
}

my @version = eval { PDL::Transform::Proj4::proj_version() };
is $@, '', 'proj_version no die';
diag "PROJ version: @version";
Expand All @@ -38,10 +29,10 @@ my $xy_exp = double [
];

my ($xy) = PDL::Transform::Proj4::fwd_transform($lonlat, $proj);
ok( tapprox( $xy, $xy_exp ) );
is_pdl $xy, $xy_exp;

my ($lonlat2) = PDL::Transform::Proj4::inv_transform($xy, $proj);
ok( tapprox( $lonlat2, $lonlat ) );
is_pdl $lonlat2, $lonlat;

# Do the corners of a cyl eq map, and see what we get...
my $cyl_eq = "+proj=eqc +lon_0=0 +datum=WGS84";
Expand All @@ -54,7 +45,7 @@ my $xy3_exp = double [
];

my ($xy3) = PDL::Transform::Proj4::fwd_transform($lonlat3, $cyl_eq);
ok( tapprox( $xy3, $xy3_exp ) );
is_pdl $xy3, $xy3_exp;

$lonlat = ((xvals( double, 35, 17 ) - 17.0) * 10.0)->cat(
(yvals( double, 35, 17 ) - 8.0) * 10.0)->mv(2,0);
Expand Down

0 comments on commit 21979bb

Please sign in to comment.