convert PDL::Fit::{LM,Linfit} to use MatrixOps not Slatec

Libtmp/Slatec/LM.pm → Libtmp/Fit/LM.pm

@@ -25,7 +25,7 @@ use warnings;
 use PDL::Core;
 use PDL::Exporter;
 use PDL::Options;
-use PDL::Slatec; # for matrix inversion
+use PDL::MatrixOps qw(lu_decomp lu_backsub inv); # for matrix inversion
 
 our @EXPORT_OK  = qw(lmfit tlmfit);
 our %EXPORT_TAGS = (Func=>\@EXPORT_OK);
@@ -107,7 +107,7 @@ sub PDL::lmfit {
   my ($maxiter,$eps) = map {$opt->{$_}} qw/Maxiter Eps/;
   # initialize some variables
   my ($isig2,$chisq) = (1/($sig*$sig),0); #$isig2="inverse of sigma squared"
-  my ($ym,$al,$cov,$bet,$oldbet,$olda,$oldal,$ochisq,$di,$pivt,$info) =
+  my ($ym,$al,$cov,$bet,$oldbet,$olda,$oldal,$ochisq,$di) =
     map {null} (0..10);
   my ($aldiag,$codiag);  # the diagonals for later updating
   # this will break broadcasting
@@ -120,9 +120,8 @@ sub PDL::lmfit {
       $cov .= $al;
       # local $PDL::debug = 1;
       $codiag .= $aldiag*(1+$lambda);
-      gefa $cov, $pivt, $info;     # gefa + gesl = solution by Gaussian elem.
-      gesl $cov, $pivt, $bet, 0;   # solution returned in $bet
-      # lusd($cov,$bet,$da);
+      my ($lu, $perm, $par) = lu_decomp($cov);
+      $bet .= lu_backsub($lu,$perm,$par, $bet);
       # print "changing by $da\n";
       $c += $bet;                  # what we used to call $da is now $bet
     }
@@ -155,7 +154,7 @@ sub PDL::lmfit {
   } while ($iter++==0 || $iter < $maxiter && $di/$chisq > $eps);
   barf "iteration did not converge" if $iter >= $maxiter && $di/$chisq > $eps;
   # return inv $al as estimate of covariance matrix
-  return wantarray ? ($ym,$c,matinv($al),$iter) : $ym;
+  return wantarray ? ($ym,$c,inv($al),$iter) : $ym;
 }
 *lmfit = \&PDL::lmfit;
 

Libtmp/Slatec/Linfit.pm → Libtmp/Fit/Linfit.pm

@@ -113,7 +113,7 @@ use PDL::Core;
 use PDL::Basic;
 use PDL::Exporter;
 use PDL::Options ':Func';
-use PDL::Slatec; # For matinv()
+use PDL::MatrixOps qw(inv);
 
 our @EXPORT_OK  = qw( linfit1d );
 our %EXPORT_TAGS = (Func=>\@EXPORT_OK);
@@ -138,7 +138,7 @@ sub PDL::linfit1d {
    my $C = $M->transpose x ($M * $wt->dummy(0)) ;
    my $Y = $M->transpose x ($y2->dummy(0) * $wt->dummy(0));
    # Fitted coefficients vector
-   my $c = matinv($C) x $Y;
+   my $c = inv($C) x $Y;
    # Fitted data
    my $yfit = ($M x $c)->clump(2) * $ymean; # Remove first dim=1, un-normalise
    return $yfit if !wantarray;

Libtmp/Fit/Makefile.PL

@@ -1,6 +1,7 @@
 use strict;
 use warnings;
 use ExtUtils::MakeMaker;
+use PDL::Core::Dev;
 
 WriteMakefile(
   NAME => 'PDL::Fit',

Libtmp/Slatec/Makefile.PL

@@ -96,8 +96,6 @@ if ($Config{osname} =~ /darwin/ && version->parse($Config{osvers}) >= version->p
        $hash{LDDLFLAGS} =~ s/-arch i386/ /g;
 }
 
-my @fit = qw(Linfit.pm LM.pm);
-@{$hash{PM}}{@fit} = map '$(INST_LIBDIR)/Fit/'.$_, @fit;
 my @filter = qw(LinPred.pm);
 @{$hash{PM}}{@filter} = map '$(INST_LIBDIR)/Filter/'.$_, @filter;
 my @gauss = qw(Gaussian.pm);

MANIFEST

@@ -269,8 +269,13 @@ Libtmp/Fit/Gaussian/gauss.c
 Libtmp/Fit/Gaussian/gaussian.pd
 Libtmp/Fit/Gaussian/Makefile.PL
 Libtmp/Fit/Gaussian/t/gauss.t
+Libtmp/Fit/Linfit.pm
+Libtmp/Fit/LM.pm
+Libtmp/Fit/lmfit_example.pl
 Libtmp/Fit/Makefile.PL
 Libtmp/Fit/Polynomial.pm
+Libtmp/Fit/t/linfit.t
+Libtmp/Fit/t/lm.t
 Libtmp/Fit/t/poly.t
 Libtmp/Func.pm
 Libtmp/Func_demo.pm
@@ -281,10 +286,7 @@ Libtmp/Simplex/Simplex.pm
 Libtmp/Simplex/t/simplex.t
 Libtmp/Slatec/barf.c
 Libtmp/Slatec/Gaussian.pm
-Libtmp/Slatec/Linfit.pm
 Libtmp/Slatec/LinPred.pm
-Libtmp/Slatec/LM.pm
-Libtmp/Slatec/lmfit_example.pl
 Libtmp/Slatec/Makefile.PL
 Libtmp/Slatec/slatec.pd
 Libtmp/Slatec/slatec/bvalu.f
@@ -400,8 +402,6 @@ Libtmp/Slatec/slatec/xermsg.f
 Libtmp/Slatec/slatec/xerprn.f
 Libtmp/Slatec/slatec/xersve.f
 Libtmp/Slatec/slatec/xgetua.f
-Libtmp/Slatec/t/linfit.t
-Libtmp/Slatec/t/lm.t
 Libtmp/Slatec/t/slatec-polyfit-weight.t
 Libtmp/Slatec/t/slatec.t
 Libtmp/t/func.t