1# -*-perl-*-
2
3use strict;
4use warnings;
5use Test::More;
6use PDL::LiteF;
7use PDL::Config;
8use PDL::IO::Misc 'rcols';
9
10BEGIN {
11    if ($PDL::Config{WITH_SLATEC}) {
12	eval " use PDL::Fit::LM; ";
13	unless ($@) {
14	    plan tests => 2;
15	}
16	else {
17	    plan skip_all => 'PDL::Fit::LM did not load. Is PDL::Slatec available?';
18	}
19    }
20    else {
21	plan skip_all => 'PDL::Fit::LM not available (needs PDL::Slatec)';
22    }
23}
24
25my ($t,$count,$sigma)=rcols(\*DATA,0,1,2);
26my $initp = pdl(10,900,80,27,225);
27my $gnuplot_pf_unweighted = pdl(7.96, 282.5, 70.0, 28.5, 117.7);
28my $gnuplot_pf_weighted = pdl(5.53, 290.7, 46.6, 33.3, 162.7);
29
30my ($yf,$pf,$cf,$if) = lmfit($t, $count, 1, \&const_2exp, $initp);
31ok(all(abs(log10($pf/$gnuplot_pf_unweighted))<0.02),"Unweighted fit");
32
33($yf,$pf,$cf,$if) = lmfit($t, $count, $sigma, \&const_2exp, $initp);
34ok(all(abs(log10($pf/$gnuplot_pf_weighted))<0.02),"Weighted fit");
35
361;
37
38sub const_2exp{
39#constant plus 2 exponentials
40    my ($x,$par,$ym,$dyda) = @_;
41    my ($a1,$a2,$a3,$a4,$a5) = map { $par->slice("($_)") } (0..4);
42    $ym .= $a1 + $a2*exp(-$x/$a4) + $a3*exp(-$x/$a5);
43    my (@dy) = map {$dyda -> slice(",($_)") } (0..4);
44    $dy[0] .= 1;
45    $dy[1] .= exp(-$x/$a4);
46    $dy[2] .= exp(-$x/$a5);
47    $dy[3] .= $a2 * $x * exp(-$x/$a4)/$a4/$a4;
48    $dy[4] .= $a3 * $x * exp(-$x/$a5)/$a5/$a5;
49}
50
51
52__DATA__
53# $Id: silver.dat,v 1.1.1.1 1998/04/15 19:16:42 lhecking Exp $
54# This sample data was distributed with Gnuplot, which contains the following notice:
55# Copyright (C) 1986 - 1993, 1998, 2004, 2007  Thomas Williams, Colin Kelley
56# Permission to use, copy, and distribute this software and its
57# documentation for any purpose with or without fee is hereby granted,
58# provided that the above copyright notice appear in all copies and that
59# both that copyright notice and this permission notice appear in
60# supporting documentation.
6110.000000 280.000000 16.733201
6220.000000 191.000000 13.820275
6330.000000 152.000000 12.328828
6440.000000 150.000000 12.247449
6550.000000 104.000000 10.198039
6660.000000 77.000000 8.774964
6770.000000 69.000000 8.306624
6880.000000 60.000000 7.745967
6990.000000 60.000000 7.745967
70100.000000 51.000000 7.141428
71110.000000 41.000000 6.403124
72120.000000 34.000000 5.830952
73130.000000 35.000000 5.916080
74140.000000 34.000000 5.830952
75150.000000 24.000000 4.898979
76160.000000 24.000000 4.898979
77170.000000 19.000000 4.358899
78180.000000 21.000000 4.582576
79190.000000 20.000000 4.472136
80200.000000 18.000000 4.242641
81210.000000 21.000000 4.582576
82220.000000 15.000000 3.872983
83230.000000 19.000000 4.358899
84240.000000 12.000000 3.464102
85250.000000 20.000000 4.472136
86260.000000 20.000000 4.472136
87270.000000 18.000000 4.242641
88280.000000 18.000000 4.242641
89290.000000 20.000000 4.472136
90300.000000 12.000000 3.464102
91310.000000 26.000000 5.099020
92320.000000 17.000000 4.123106
93330.000000 8.000000 2.828427
94340.000000 6.000000 2.449490
95350.000000 8.000000 2.828427
96360.000000 10.000000 3.162278
97370.000000 20.000000 4.472136
98380.000000 14.000000 3.741657
99390.000000 8.000000 2.828427
100400.000000 10.000000 3.162278
101410.000000 9.000000 3.000000
102420.000000 8.000000 2.828427
103430.000000 10.000000 3.162278
104440.000000 13.000000 3.605551
105450.000000 9.000000 3.000000
106460.000000 5.000000 2.236068
107470.000000 7.000000 2.645751
108480.000000 11.000000 3.316625
109500.000000 7.000000 2.645751
110510.000000 9.000000 3.000000
111520.000000 12.000000 3.464102
112530.000000 4.000000 2.000000
113540.000000 7.000000 2.645751
114550.000000 10.000000 3.162278
115560.000000 9.000000 3.000000
116580.000000 8.000000 2.828427
117590.000000 9.000000 3.000000
118600.000000 5.000000 2.236068
119