LM.pm (PDL-2.074) | : | LM.pm (PDL-2.075) | ||
---|---|---|---|---|

skipping to change at line 113 | skipping to change at line 113 | |||

sub PDL::lmfit { | sub PDL::lmfit { | |||

my ($x,$y,$sig,$func,$c,$opt) = @_; # not using $ia right now | my ($x,$y,$sig,$func,$c,$opt) = @_; # not using $ia right now | |||

$opt = {iparse( { Maxiter => 200, | $opt = {iparse( { Maxiter => 200, | |||

Eps => 1e-4}, ifhref($opt))}; | Eps => 1e-4}, ifhref($opt))}; | |||

my ($maxiter,$eps) = map {$opt->{$_}} qw/Maxiter Eps/; | my ($maxiter,$eps) = map {$opt->{$_}} qw/Maxiter Eps/; | |||

# initialize some variables | # initialize some variables | |||

my ($isig2,$chisq) = (1/($sig*$sig),0); #$isig2="inverse of sigma squared" | 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,$pivt,$info) = | |||

map {null} (0..10); | map {null} (0..10); | |||

my ($aldiag,$codiag); # the diagonals for later updating | my ($aldiag,$codiag); # the diagonals for later updating | |||

# this will break threading | # this will break broadcasting | |||

my $dyda = zeroes($x->type,$x->getdim(0),$c->getdim(0)); | my $dyda = zeroes($x->type,$x->getdim(0),$c->getdim(0)); | |||

my $alv = zeroes($x->type,$x->getdim(0),$c->getdim(0),$c->getdim(0)); | my $alv = zeroes($x->type,$x->getdim(0),$c->getdim(0),$c->getdim(0)); | |||

my ($iter,$lambda) = (0,0.001); | my ($iter,$lambda) = (0,0.001); | |||

do { | do { | |||

if ($iter>0) { | if ($iter>0) { | |||

$cov .= $al; | $cov .= $al; | |||

# local $PDL::debug = 1; | # local $PDL::debug = 1; | |||

$codiag .= $aldiag*(1+$lambda); | $codiag .= $aldiag*(1+$lambda); | |||

gefa $cov, $pivt, $info; # gefa + gesl = solution by Gaussian elem. | gefa $cov, $pivt, $info; # gefa + gesl = solution by Gaussian elem. | |||

skipping to change at line 262 | skipping to change at line 262 | |||

} | } | |||

=cut | =cut | |||

# the OtherPar is the sub routine ref | # the OtherPar is the sub routine ref | |||

=head2 tlmfit | =head2 tlmfit | |||

=for ref | =for ref | |||

threaded version of Levenberg-Marquardt fitting routine mfit | broadcasted version of Levenberg-Marquardt fitting routine mfit | |||

=for example | =for example | |||

tlmfit $x, $y, float(1)->dummy(0), $na, float(200), float(1e-4), | tlmfit $x, $y, float(1)->dummy(0), $na, float(200), float(1e-4), | |||

$ym=null, $afit=null, \&expdec; | $ym=null, $afit=null, \&expdec; | |||

=for sig | =for sig | |||

Signature: tlmfit(x(n);y(n);sigma(n);initp(m);iter();eps();[o] ym(n);[o] final p(m); | Signature: tlmfit(x(n);y(n);sigma(n);initp(m);iter();eps();[o] ym(n);[o] final p(m); | |||

OtherPar => subref) | OtherPar => subref) | |||

a threaded version of C<lmfit> by using perl threading. Direct | a broadcasted version of C<lmfit> by using perl broadcasting. Direct | |||

threading in C<lmfit> seemed difficult since we have an if condition | broadcasting in C<lmfit> seemed difficult since we have an if condition | |||

in the iteration. In principle that can be worked around by | in the iteration. In principle that can be worked around by | |||

using C<where> but .... Send a threaded C<lmfit> version if | using C<where> but .... Send a broadcasted C<lmfit> version if | |||

you work it out! | you work it out! | |||

Since we are using perl threading here speed is not really great | Since we are using perl broadcasting here speed is not really great | |||

but it is just convenient to have a threaded version for many | but it is just convenient to have a broadcasted version for many | |||

applications (no explicit for-loops required, etc). Suffers from | applications (no explicit for-loops required, etc). Suffers from | |||

some of the current limitations of perl level threading. | some of the current limitations of perl level broadcasting. | |||

=cut | =cut | |||

thread_define 'tlmfit(x(n);y(n);sigma(n);initp(m);iter();eps();[o] ym(n);[o] fin alp(m)), | broadcast_define 'tlmfit(x(n);y(n);sigma(n);initp(m);iter();eps();[o] ym(n);[o] finalp(m)), | |||

NOtherPars => 1', | NOtherPars => 1', | |||

over { | over { | |||

$_[7] .= $_[3]; # copy our parameter guess into the output | $_[7] .= $_[3]; # copy our parameter guess into the output | |||

$_[6] .= PDL::lmfit $_[0],$_[1],$_[2],$_[8],$_[7],{Maxiter => $_[4], | $_[6] .= PDL::lmfit $_[0],$_[1],$_[2],$_[8],$_[7],{Maxiter => $_[4], | |||

Eps => $_[5]}; | Eps => $_[5]}; | |||

}; | }; | |||

1; | 1; | |||

=head1 BUGS | =head1 BUGS | |||

End of changes. 7 change blocks. | ||||

9 lines changed or deleted | | 9 lines changed or added |