t/Multifit.t

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