lib/Math/GSL/Multifit/Test.pm

   1 package Math::GSL::Multifit::Test;
   2 use base q{Test::Class};
   3 use Test::More;
   4 use Math::GSL::Multifit qw/:all/;
   5 use Math::GSL::Matrix qw/:all/;
   6 use Math::GSL::Vector qw/:all/;
   7 use Math::GSL::BLAS qw/:all/;
   8 use Math::GSL::Machine qw/:all/;
   9 use Math::GSL::Errno qw/:all/;
  10 use Math::GSL qw/:all/;
  11 use Data::Dumper;
  12 use strict;
  13
  14 sub make_fixture : Test(setup) {
  15 }
  16
  17 sub teardown : Test(teardown) {
  18 }
  19
  20 sub GSL_MULTIFIT_LINEAR_ALLOC : Tests {
  21  my $multi = gsl_multifit_linear_alloc(2,2);
  22  isa_ok($multi, 'Math::GSL::Multifit');
  23 }
  24
  25 sub GSL_MULTIFIT_LINEAR_EST : Tests {
  26  my $M = Math::GSL::Matrix->new(5,5);
  27  $M->set_row(0, [4.271520, -0.526675,  0.957930,  0.267750, -0.103610]);
  28  $M->set_row(1, [-0.526675,  5.701680, -0.098080,  0.641845,  0.429780]);
  29  $M->set_row(2, [0.957930, -0.098080,  4.584790,  0.375865,  1.510810]);
  30  $M->set_row(3, [0.267750,  0.641845,  0.375865,  4.422720,  0.392210]);
  31  $M->set_row(4, [-0.103610,  0.429780,  1.510810,  0.392210,  5.782750]);
  32  my $c = Math::GSL::Vector->new([-0.627020,   0.848674,   0.216877,  -0.057883,   0.596668]);
  33  my $x = Math::GSL::Vector->new([0.99932,   0.23858,   0.19797,   1.44008,  -0.15335]);
  34  my @got = gsl_multifit_linear_est($x->raw, $c->raw, $M->raw);
  35  is($got[0],0);
  36  ok(is_similar_relative($got[1], -5.56037032230000*(10**-1), 256*$GSL_DBL_EPSILON));
  37  ok(is_similar_relative($got[2], 3.91891123349318, 256*$GSL_DBL_EPSILON));
  38 }
  39
  40 sub GSL_MULTIFIT_LINEAR__GSL_MULTIFIT_LINEAR_RESIDUAL : Tests {
  41
  42 my @filip_x = ( -6.860120914, -4.324130045, -4.358625055,
  43 -4.358426747, -6.955852379, -6.661145254, -6.355462942, -6.118102026,
  44 -7.115148017, -6.815308569, -6.519993057, -6.204119983, -5.853871964,
  45 -6.109523091, -5.79832982, -5.482672118, -5.171791386, -4.851705903,
  46 -4.517126416, -4.143573228, -3.709075441, -3.499489089, -6.300769497,
  47 -5.953504836, -5.642065153, -5.031376979, -4.680685696, -4.329846955,
  48 -3.928486195, -8.56735134, -8.363211311, -8.107682739, -7.823908741,
  49 -7.522878745, -7.218819279, -6.920818754, -6.628932138, -6.323946875,
  50 -5.991399828, -8.781464495, -8.663140179, -8.473531488, -8.247337057,
  51 -7.971428747, -7.676129393, -7.352812702, -7.072065318, -6.774174009,
  52 -6.478861916, -6.159517513, -6.835647144, -6.53165267, -6.224098421,
  53 -5.910094889, -5.598599459, -5.290645224, -4.974284616, -4.64454848,
  54 -4.290560426, -3.885055584, -3.408378962, -3.13200249, -8.726767166,
  55 -8.66695597, -8.511026475, -8.165388579, -7.886056648, -7.588043762,
  56 -7.283412422, -6.995678626, -6.691862621, -6.392544977, -6.067374056,
  57 -6.684029655, -6.378719832, -6.065855188, -5.752272167, -5.132414673,
  58 -4.811352704, -4.098269308, -3.66174277, -3.2644011);
  59
  60 my $filip_y = [0.8116, 0.9072, 0.9052, 0.9039, 0.8053, 0.8377,
  61 0.8667, 0.8809, 0.7975, 0.8162, 0.8515, 0.8766, 0.8885, 0.8859,
  62 0.8959, 0.8913, 0.8959, 0.8971, 0.9021, 0.909, 0.9139, 0.9199, 0.8692,
  63 0.8872, 0.89, 0.891, 0.8977, 0.9035, 0.9078, 0.7675, 0.7705, 0.7713,
  64 0.7736, 0.7775, 0.7841, 0.7971, 0.8329, 0.8641, 0.8804, 0.7668,
  65 0.7633, 0.7678, 0.7697, 0.77, 0.7749, 0.7796, 0.7897, 0.8131, 0.8498,
  66 0.8741, 0.8061, 0.846, 0.8751, 0.8856, 0.8919, 0.8934, 0.894, 0.8957,
  67 0.9047, 0.9129, 0.9209, 0.9219, 0.7739, 0.7681, 0.7665, 0.7703,
  68 0.7702, 0.7761, 0.7809, 0.7961, 0.8253, 0.8602, 0.8809, 0.8301,
  69 0.8664, 0.8834, 0.8898, 0.8964, 0.8963, 0.9074, 0.9119, 0.9228] ;
  70         my $filip_n = 82;
  71         my $filip_p = 11;
  72    my $work = gsl_multifit_linear_alloc ($filip_n, $filip_p);
  73
  74     my $X = gsl_matrix_alloc ($filip_n, $filip_p);
  75     my $y = Math::GSL::Vector->new($filip_y);
  76     my $c = gsl_vector_alloc ($filip_p);
  77     my $cov = gsl_matrix_alloc ($filip_p, $filip_p);
  78     my $r = gsl_vector_alloc($filip_n);
  79     my @expected_c = ( -1467.48961422980,
  80                               -2772.17959193342,
  81                               -2316.37108160893,
  82                               -1127.97394098372,
  83                               -354.478233703349,
  84                               -75.1242017393757,
  85                               -10.8753180355343,
  86                               -1.06221498588947,
  87                               -0.670191154593408*(10**-1),
  88                               -0.246781078275479*(10**-2),
  89                               -0.402962525080404*(10**-4) );
  90
  91     my @expected_sd  = ( 298.084530995537,
  92                                559.779865474950,
  93                                466.477572127796,
  94                                227.204274477751,
  95                                71.6478660875927,
  96                                15.2897178747400,
  97                                2.23691159816033,
  98                                0.221624321934227,
  99                                0.142363763154724*(10**-1),
 100                                0.535617408889821*(10**-3),
 101                                0.896632837373868*(10**-5) );
 102
 103     my $expected_chisq = 0.795851382172941*(10**-3);
 104     my ($i, $j);
 105     for ($i = 0 ; $i < $filip_n; $i++)
 106       {
 107         for ($j = 0; $j < $filip_p; $j++)
 108           {
 109             gsl_matrix_set($X, $i, $j, ($filip_x[$i])**$j);
 110           }
 111       }
 112
 113     my @got = gsl_multifit_linear($X, $y->raw, $c, $cov, $work);
 114     is ($got[0], 0);
 115     ok(is_similar_relative($got[1], $expected_chisq, 10**-7), "filip gsl_fit_multilinear chisq");
 116     ok(is_similar_relative(gsl_vector_get($c,0), $expected_c[0], 10**-7), "filip gsl_fit_multilinear c0") ;
 117     ok(is_similar_relative(gsl_vector_get($c,1), $expected_c[1], 10**-7), "filip gsl_fit_multilinear c1") ;
 118     ok(is_similar_relative(gsl_vector_get($c,2), $expected_c[2], 10**-7), "filip gsl_fit_multilinear c2") ;
 119     ok(is_similar_relative(gsl_vector_get($c,3), $expected_c[3], 10**-7), "filip gsl_fit_multilinear c3") ;
 120     ok(is_similar_relative(gsl_vector_get($c,4), $expected_c[4], 10**-7), "filip gsl_fit_multilinear c4") ;
 121     ok(is_similar_relative(gsl_vector_get($c,5), $expected_c[5], 10**-7), "filip gsl_fit_multilinear c5") ;
 122     ok(is_similar_relative(gsl_vector_get($c,6), $expected_c[6], 10**-7), "filip gsl_fit_multilinear c6") ;
 123     ok(is_similar_relative(gsl_vector_get($c,7), $expected_c[7], 10**-7), "filip gsl_fit_multilinear c7") ;
 124     ok(is_similar_relative(gsl_vector_get($c,8), $expected_c[8], 10**-7), "filip gsl_fit_multilinear c8") ;
 125     ok(is_similar_relative(gsl_vector_get($c,9), $expected_c[9], 10**-7), "filip gsl_fit_multilinear c9") ;
 126     ok(is_similar_relative(gsl_vector_get($c,10), $expected_c[10], 10**-7), "filip gsl_fit_multilinear c10") ;
 127
 128     my $diag = gsl_matrix_diagonal($cov);
 129
 130     ok(is_similar_relative(gsl_vector_get($diag->{vector},0), $expected_sd[0]**2, 10**-6), "filip gsl_fit_multilinear cov00") ;
 131     ok(is_similar_relative(gsl_vector_get($diag->{vector},1), $expected_sd[1]**2, 10**-6), "filip gsl_fit_multilinear cov11") ;
 132     ok(is_similar_relative(gsl_vector_get($diag->{vector},2), $expected_sd[2]**2, 10**-6), "filip gsl_fit_multilinear cov22") ;
 133     ok(is_similar_relative(gsl_vector_get($diag->{vector},3), $expected_sd[3]**2, 10**-6), "filip gsl_fit_multilinear cov33") ;
 134     ok(is_similar_relative(gsl_vector_get($diag->{vector},4), $expected_sd[4]**2, 10**-6), "filip gsl_fit_multilinear cov44") ;
 135     ok(is_similar_relative(gsl_vector_get($diag->{vector},5), $expected_sd[5]**2, 10**-6), "filip gsl_fit_multilinear cov55") ;
 136     ok(is_similar_relative(gsl_vector_get($diag->{vector},6), $expected_sd[6]**2, 10**-6), "filip gsl_fit_multilinear cov66") ;
 137     ok(is_similar_relative(gsl_vector_get($diag->{vector},7), $expected_sd[7]**2, 10**-6), "filip gsl_fit_multilinear cov77") ;
 138     ok(is_similar_relative(gsl_vector_get($diag->{vector},8), $expected_sd[8]**2, 10**-6), "filip gsl_fit_multilinear cov88") ;
 139     ok(is_similar_relative(gsl_vector_get($diag->{vector},9), $expected_sd[9]**2, 10**-6), "filip gsl_fit_multilinear cov99") ;
 140     ok(is_similar_relative(gsl_vector_get($diag->{vector},10), $expected_sd[10]**2, 10**-6), "filip gsl_fit_multilinear cov1010") ;
 141
 142     gsl_multifit_linear_residuals($X, $y->raw, $c, $r);
 143     @got = gsl_blas_ddot($r, $r);
 144     ok(is_similar_relative($got[1], $expected_chisq, 10**-7), "filip gsl_fit_multilinear residuals") ;
 145 }
 146
 147 sub GSL_MULTIFIT_WLINEAR : Tests {
 148     my $filip_n = 82;
 149     my $filip_p = 11;
 150 my @filip_x = ( -6.860120914, -4.324130045, -4.358625055,
 151 -4.358426747, -6.955852379, -6.661145254, -6.355462942, -6.118102026,
 152 -7.115148017, -6.815308569, -6.519993057, -6.204119983, -5.853871964,
 153 -6.109523091, -5.79832982, -5.482672118, -5.171791386, -4.851705903,
 154 -4.517126416, -4.143573228, -3.709075441, -3.499489089, -6.300769497,
 155 -5.953504836, -5.642065153, -5.031376979, -4.680685696, -4.329846955,
 156 -3.928486195, -8.56735134, -8.363211311, -8.107682739, -7.823908741,
 157 -7.522878745, -7.218819279, -6.920818754, -6.628932138, -6.323946875,
 158 -5.991399828, -8.781464495, -8.663140179, -8.473531488, -8.247337057,
 159 -7.971428747, -7.676129393, -7.352812702, -7.072065318, -6.774174009,
 160 -6.478861916, -6.159517513, -6.835647144, -6.53165267, -6.224098421,
 161 -5.910094889, -5.598599459, -5.290645224, -4.974284616, -4.64454848,
 162 -4.290560426, -3.885055584, -3.408378962, -3.13200249, -8.726767166,
 163 -8.66695597, -8.511026475, -8.165388579, -7.886056648, -7.588043762,
 164 -7.283412422, -6.995678626, -6.691862621, -6.392544977, -6.067374056,
 165 -6.684029655, -6.378719832, -6.065855188, -5.752272167, -5.132414673,
 166 -4.811352704, -4.098269308, -3.66174277, -3.2644011);
 167
 168 my $filip_y = [0.8116, 0.9072, 0.9052, 0.9039, 0.8053, 0.8377,
 169 0.8667, 0.8809, 0.7975, 0.8162, 0.8515, 0.8766, 0.8885, 0.8859,
 170 0.8959, 0.8913, 0.8959, 0.8971, 0.9021, 0.909, 0.9139, 0.9199, 0.8692,
 171 0.8872, 0.89, 0.891, 0.8977, 0.9035, 0.9078, 0.7675, 0.7705, 0.7713,
 172 0.7736, 0.7775, 0.7841, 0.7971, 0.8329, 0.8641, 0.8804, 0.7668,
 173 0.7633, 0.7678, 0.7697, 0.77, 0.7749, 0.7796, 0.7897, 0.8131, 0.8498,
 174 0.8741, 0.8061, 0.846, 0.8751, 0.8856, 0.8919, 0.8934, 0.894, 0.8957,
 175 0.9047, 0.9129, 0.9209, 0.9219, 0.7739, 0.7681, 0.7665, 0.7703,
 176 0.7702, 0.7761, 0.7809, 0.7961, 0.8253, 0.8602, 0.8809, 0.8301,
 177 0.8664, 0.8834, 0.8898, 0.8964, 0.8963, 0.9074, 0.9119, 0.9228] ;
 178     my $work = gsl_multifit_linear_alloc ($filip_n, $filip_p);
 179     my $X = gsl_matrix_alloc ($filip_n, $filip_p);
 180     my $y = Math::GSL::Vector->new($filip_y);
 181     my $w = gsl_vector_alloc($filip_n);
 182     my $c = gsl_vector_alloc($filip_p);
 183     my $r = gsl_vector_alloc($filip_n);
 184     my $cov = gsl_matrix_alloc($filip_p, $filip_p);
 185
 186     my @expected_c = ( -1467.48961422980,
 187                               -2772.17959193342,
 188                               -2316.37108160893,
 189                               -1127.97394098372,
 190                               -354.478233703349,
 191                               -75.1242017393757,
 192                               -10.8753180355343,
 193                               -1.06221498588947,
 194                               -0.670191154593408*(10**-1),
 195                               -0.246781078275479*(10**-2),
 196                               -0.402962525080404*(10**-4) );
 197
 198     # computed using GNU Calc
 199
 200     my @expected_cov =( [  7.9269341767252183262588583867942e9,  1.4880416622254098343441063389706e10, 1.2385811858111487905481427591107e10, 6.0210784406215266653697715794241e9, 1.8936652526181982747116667336389e9, 4.0274900618493109653998118587093e8, 5.8685468011819735806180092394606e7, 5.7873451475721689084330083708901e6,  3.6982719848703747920663262917032e5,  1.3834818802741350637527054170891e4,   2.301758578713219280719633494302e2  ],
 201       [ 1.4880416622254098334697515488559e10, 2.7955091668548290835529555438088e10, 2.3286604504243362691678565997033e10, 1.132895006796272983689297219686e10, 3.5657281653312473123348357644683e9, 7.5893300392314445528176646366087e8, 1.1066654886143524811964131660002e8, 1.0921285448484575110763947787775e7,  6.9838139975394769253353547606971e5,  2.6143091775349597218939272614126e4,  4.3523386330348588614289505633539e2  ],
 202       [ 1.2385811858111487890788272968677e10, 2.3286604504243362677757802422747e10, 1.9412787917766676553608636489674e10, 9.4516246492862131849077729250098e9, 2.9771226694709917550143152097252e9, 6.3413035086730038062129508949859e8, 9.2536164488309401636559552742339e7, 9.1386304643423333815338760248027e6,  5.8479478338916429826337004060941e5,  2.1905933113294737443808429764554e4,  3.6493161325305557266196635180155e2  ],
 203       [ 6.0210784406215266545770691532365e9,  1.1328950067962729823273441573365e10, 9.4516246492862131792040001429636e9,  4.6053152992000107509329772255094e9, 1.4517147860312147098138030287038e9, 3.0944988323328589376402579060072e8, 4.5190223822292688669369522708712e7, 4.4660958693678497534529855690752e6,  2.8599340736122198213681258676423e5,  1.0720394998549386596165641244705e4,  1.7870937745661967319298031044424e2  ],
 204       [ 1.8936652526181982701620450132636e9,  3.5657281653312473058825073094524e9,  2.9771226694709917514149924058297e9,  1.451714786031214708936087401632e9,  4.5796563896564815123074920050827e8, 9.7693972414561515534525103622773e7, 1.427717861635658545863942948444e7,  1.4120161287735817621354292900338e6,  9.0484361228623960006818614875557e4,   3.394106783764852373199087455398e3,  5.6617406468519495376287407526295e1  ],
 205     [ 4.0274900618493109532650887473599e8,   7.589330039231444534478894935778e8,  6.3413035086730037947153564986653e8,   3.09449883233285893390542947998e8,  9.7693972414561515475770399055121e7, 2.0855726248311948992114244257719e7, 3.0501263034740400533872858749566e6, 3.0187475839310308153394428784224e5,  1.9358204633534233524477930175632e4,  7.2662989867560017077361942813911e2,  1.2129002231061036467607394277965e1  ],
 206       [  5.868546801181973559370854830868e7,  1.1066654886143524778548044386795e8,  9.2536164488309401413296494869777e7,  4.5190223822292688587853853162072e7, 1.4277178616356585441556046753562e7, 3.050126303474040051574715592746e6,  4.4639982579046340884744460329946e5, 4.4212093985989836047285007760238e4,  2.8371395028774486687625333589972e3,  1.0656694507620102300567296504381e2,  1.7799982046359973175080475654123e0  ],
 207       [ 5.7873451475721688839974153925406e6,  1.0921285448484575071271480643397e7,  9.1386304643423333540728480344578e6,  4.4660958693678497427674903565664e6, 1.4120161287735817596182229182587e6, 3.0187475839310308117812257613082e5, 4.4212093985989836021482392757677e4, 4.3818874017028389517560906916315e3,   2.813828775753142855163154605027e2,  1.0576188138416671883232607188969e1,  1.7676976288918295012452853715408e-1 ],
 208       [ 3.6982719848703747742568351456818e5,  6.9838139975394768959780068745979e5,  5.8479478338916429616547638954781e5,  2.8599340736122198128717796825489e5, 9.0484361228623959793493985226792e4, 1.9358204633534233490579641064343e4, 2.8371395028774486654873647731797e3, 2.8138287757531428535592907878017e2,  1.8081118503579798222896804627964e1,  6.8005074291434681866415478598732e-1, 1.1373581557749643543869665860719e-2 ],
 209       [ 1.3834818802741350562839757244708e4,   2.614309177534959709397445440919e4,  2.1905933113294737352721470167247e4,  1.0720394998549386558251721913182e4, 3.3941067837648523632905604575131e3, 7.2662989867560016909534954790835e2, 1.0656694507620102282337905013451e2, 1.0576188138416671871337685672492e1,  6.8005074291434681828743281967838e-1, 2.5593857187900736057022477529078e-2, 4.2831487599116264442963102045936e-4 ],
 210       [ 2.3017585787132192669801658674163e2,  4.3523386330348588381716460685124e2,  3.6493161325305557094116270974735e2,  1.7870937745661967246233792737255e2, 5.6617406468519495180024059284629e1, 1.2129002231061036433003571679329e1, 1.7799982046359973135014027410646e0, 1.7676976288918294983059118597214e-1, 1.137358155774964353146460100337e-2,  4.283148759911626442000316269063e-4,  7.172253875245080423800933453952e-6  ] );
 211
 212     my $expected_chisq = 0.795851382172941*(10**-3);
 213     my ($i, $j);
 214     for ($i = 0 ; $i < $filip_n; $i++)
 215       {
 216         for ($j = 0; $j < $filip_p; $j++)
 217           {
 218             gsl_matrix_set($X, $i, $j, $filip_x[$i]**$j);
 219           }
 220       }
 221
 222     gsl_vector_set_all($w, 1.0);
 223
 224     my @got = gsl_multifit_wlinear($X, $w, $y->raw, $c, $cov, $work);
 225     ok(is_similar_relative(gsl_vector_get($c,0), $expected_c[0], 10**-7), "filip gsl_fit_multilinear c0");
 226     ok(is_similar_relative(gsl_vector_get($c,1), $expected_c[1], 10**-7), "filip gsl_fit_multilinear c1") ;
 227     ok(is_similar_relative(gsl_vector_get($c,2), $expected_c[2], 10**-7), "filip gsl_fit_multilinear c2") ;
 228     ok(is_similar_relative(gsl_vector_get($c,3), $expected_c[3], 10**-7), "filip gsl_fit_multilinear c3") ;
 229     ok(is_similar_relative(gsl_vector_get($c,4), $expected_c[4], 10**-7), "filip gsl_fit_multilinear c4") ;
 230     ok(is_similar_relative(gsl_vector_get($c,5), $expected_c[5], 10**-7), "filip gsl_fit_multilinear c5") ;
 231     ok(is_similar_relative(gsl_vector_get($c,6), $expected_c[6], 10**-7), "filip gsl_fit_multilinear c6") ;
 232     ok(is_similar_relative(gsl_vector_get($c,7), $expected_c[7], 10**-7), "filip gsl_fit_multilinear c7") ;
 233     ok(is_similar_relative(gsl_vector_get($c,8), $expected_c[8], 10**-7), "filip gsl_fit_multilinear c8") ;
 234     ok(is_similar_relative(gsl_vector_get($c,9), $expected_c[9], 10**-7), "filip gsl_fit_multilinear c9") ;
 235     ok(is_similar_relative(gsl_vector_get($c,10), $expected_c[10], 10**-7), "filip gsl_fit_multilinear c10") ;
 236
 237
 238     for ($i = 0; $i < $filip_p; $i++)
 239       {
 240         for ($j = 0; $j < $filip_p; $j++)
 241           {
 242             ok(is_similar_relative(gsl_matrix_get($cov,$i,$j), $expected_cov[$i][$j], 10**-6),"filip gsl_fit_wmultilinear cov($i,$j)");
 243           }
 244       }
 245
 246     ok(is_similar_relative($got[1], $expected_chisq, 10**-7), "filip gsl_fit_wmultilinear chisq") ;
 247
 248     gsl_multifit_linear_residuals($X, $y->raw, $c, $r);
 249     @got = gsl_blas_ddot($r, $r);
 250     ok(is_similar_relative($got[1], $expected_chisq, 10**-7), "filip gsl_fit_wmultilinear residuals") ;
 251 }
 252
 253 sub GSL_MULTIFIT_FSOLVER_ALLOC : Tests {
 254    local $TODO = "the only gsl_multifit_fsolver_type is commented out in gsl_multifit_nlin.h";
 255    ok(0);
 256    #my $fsolver = gsl_multifit_fsolver_alloc($gsl_multifit_fsolver_gradient, 100, 3);
 257    #isa_ok($fsolver, 'Math::GSL::Multifit');
 258 }
 259 1;