htswanalysis/MACS/lib/gsl/gsl-1.11/multifit/test_hahn1.c

   1 const size_t hahn1_N = 236;
   2 const size_t hahn1_P = 7;
   3
   4 /* double hahn1_x0[7] = { 10, -1, 0.05, -0.00001, -0.05, 0.001, -0.000001 }; */
   5
   6 double hahn1_x0[7] = { 1, -0.1, 0.005, -0.000001, -0.005, 0.0001, -0.0000001};
   7
   8 double hahn1_x[7] = {
   9 1.0776351733E+00,
  10 -1.2269296921E-01,
  11 4.0863750610E-03,
  12 -1.4262662514E-06,
  13 -5.7609940901E-03,
  14 2.4053735503E-04,
  15 -1.2314450199E-07
  16 };
  17
  18 double hahn1_sumsq = 1.5324382854E+00;
  19
  20 double hahn1_sigma[7] = {
  21   1.7070154742E-01,
  22   1.2000289189E-02,
  23   2.2508314937E-04,
  24   2.7578037666E-07,
  25   2.4712888219E-04,
  26   1.0449373768E-05,
  27   1.3027335327E-08
  28 };
  29
  30 double hahn1_F1[236] = {
  31         .591E0,
  32        1.547E0,
  33        2.902E0,
  34        2.894E0,
  35        4.703E0,
  36        6.307E0,
  37        7.03E0 ,
  38        7.898E0,
  39        9.470E0,
  40        9.484E0,
  41       10.072E0,
  42       10.163E0,
  43       11.615E0,
  44       12.005E0,
  45       12.478E0,
  46       12.982E0,
  47       12.970E0,
  48       13.926E0,
  49       14.452E0,
  50       14.404E0,
  51       15.190E0,
  52       15.550E0,
  53       15.528E0,
  54       15.499E0,
  55       16.131E0,
  56       16.438E0,
  57       16.387E0,
  58       16.549E0,
  59       16.872E0,
  60       16.830E0,
  61       16.926E0,
  62       16.907E0,
  63       16.966E0,
  64       17.060E0,
  65       17.122E0,
  66       17.311E0,
  67       17.355E0,
  68       17.668E0,
  69       17.767E0,
  70       17.803E0,
  71       17.765E0,
  72       17.768E0,
  73       17.736E0,
  74       17.858E0,
  75       17.877E0,
  76       17.912E0,
  77       18.046E0,
  78       18.085E0,
  79       18.291E0,
  80       18.357E0,
  81       18.426E0,
  82       18.584E0,
  83       18.610E0,
  84       18.870E0,
  85       18.795E0,
  86       19.111E0,
  87         .367E0,
  88         .796E0,
  89        0.892E0,
  90        1.903E0,
  91        2.150E0,
  92        3.697E0,
  93        5.870E0,
  94        6.421E0,
  95        7.422E0,
  96        9.944E0,
  97       11.023E0,
  98       11.87E0 ,
  99       12.786E0,
 100       14.067E0,
 101       13.974E0,
 102       14.462E0,
 103       14.464E0,
 104       15.381E0,
 105       15.483E0,
 106       15.59E0 ,
 107       16.075E0,
 108       16.347E0,
 109       16.181E0,
 110       16.915E0,
 111       17.003E0,
 112       16.978E0,
 113       17.756E0,
 114       17.808E0,
 115       17.868E0,
 116       18.481E0,
 117       18.486E0,
 118       19.090E0,
 119       16.062E0,
 120       16.337E0,
 121       16.345E0,
 122       16.388E0,
 123       17.159E0,
 124       17.116E0,
 125       17.164E0,
 126       17.123E0,
 127       17.979E0,
 128       17.974E0,
 129       18.007E0,
 130       17.993E0,
 131       18.523E0,
 132       18.669E0,
 133       18.617E0,
 134       19.371E0,
 135       19.330E0,
 136        0.080E0,
 137        0.248E0,
 138        1.089E0,
 139        1.418E0,
 140        2.278E0,
 141        3.624E0,
 142        4.574E0,
 143        5.556E0,
 144        7.267E0,
 145        7.695E0,
 146        9.136E0,
 147        9.959E0,
 148        9.957E0,
 149       11.600E0,
 150       13.138E0,
 151       13.564E0,
 152       13.871E0,
 153       13.994E0,
 154       14.947E0,
 155       15.473E0,
 156       15.379E0,
 157       15.455E0,
 158       15.908E0,
 159       16.114E0,
 160       17.071E0,
 161       17.135E0,
 162       17.282E0,
 163       17.368E0,
 164       17.483E0,
 165       17.764E0,
 166       18.185E0,
 167       18.271E0,
 168       18.236E0,
 169       18.237E0,
 170       18.523E0,
 171       18.627E0,
 172       18.665E0,
 173       19.086E0,
 174        0.214E0,
 175        0.943E0,
 176        1.429E0,
 177        2.241E0,
 178        2.951E0,
 179        3.782E0,
 180        4.757E0,
 181        5.602E0,
 182        7.169E0,
 183        8.920E0,
 184       10.055E0,
 185       12.035E0,
 186       12.861E0,
 187       13.436E0,
 188       14.167E0,
 189       14.755E0,
 190       15.168E0,
 191       15.651E0,
 192       15.746E0,
 193       16.216E0,
 194       16.445E0,
 195       16.965E0,
 196       17.121E0,
 197       17.206E0,
 198       17.250E0,
 199       17.339E0,
 200       17.793E0,
 201       18.123E0,
 202       18.49E0 ,
 203       18.566E0,
 204       18.645E0,
 205       18.706E0,
 206       18.924E0,
 207       19.1E0  ,
 208        0.375E0,
 209        0.471E0,
 210        1.504E0,
 211        2.204E0,
 212        2.813E0,
 213        4.765E0,
 214        9.835E0,
 215       10.040E0,
 216       11.946E0,
 217       12.596E0,
 218       13.303E0,
 219       13.922E0,
 220       14.440E0,
 221       14.951E0,
 222       15.627E0,
 223       15.639E0,
 224       15.814E0,
 225       16.315E0,
 226       16.334E0,
 227       16.430E0,
 228       16.423E0,
 229       17.024E0,
 230       17.009E0,
 231       17.165E0,
 232       17.134E0,
 233       17.349E0,
 234       17.576E0,
 235       17.848E0,
 236       18.090E0,
 237       18.276E0,
 238       18.404E0,
 239       18.519E0,
 240       19.133E0,
 241       19.074E0,
 242       19.239E0,
 243       19.280E0,
 244       19.101E0,
 245       19.398E0,
 246       19.252E0,
 247       19.89E0 ,
 248       20.007E0,
 249       19.929E0,
 250       19.268E0,
 251       19.324E0,
 252       20.049E0,
 253       20.107E0,
 254       20.062E0,
 255       20.065E0,
 256       19.286E0,
 257       19.972E0,
 258       20.088E0,
 259       20.743E0,
 260       20.83E0 ,
 261       20.935E0,
 262       21.035E0,
 263       20.93E0 ,
 264       21.074E0,
 265       21.085E0,
 266       20.935E0
 267 };
 268
 269
 270 double hahn1_F0[236] = {
 271     24.41E0,
 272     34.82E0,
 273     44.09E0,
 274     45.07E0,
 275     54.98E0,
 276     65.51E0,
 277     70.53E0,
 278     75.70E0,
 279     89.57E0,
 280     91.14E0,
 281     96.40E0,
 282     97.19E0,
 283    114.26E0,
 284    120.25E0,
 285    127.08E0,
 286    133.55E0,
 287    133.61E0,
 288    158.67E0,
 289    172.74E0,
 290    171.31E0,
 291    202.14E0,
 292    220.55E0,
 293    221.05E0,
 294    221.39E0,
 295    250.99E0,
 296    268.99E0,
 297    271.80E0,
 298    271.97E0,
 299    321.31E0,
 300    321.69E0,
 301    330.14E0,
 302    333.03E0,
 303    333.47E0,
 304    340.77E0,
 305    345.65E0,
 306    373.11E0,
 307    373.79E0,
 308    411.82E0,
 309    419.51E0,
 310    421.59E0,
 311    422.02E0,
 312    422.47E0,
 313    422.61E0,
 314    441.75E0,
 315    447.41E0,
 316    448.7E0 ,
 317    472.89E0,
 318    476.69E0,
 319    522.47E0,
 320    522.62E0,
 321    524.43E0,
 322    546.75E0,
 323    549.53E0,
 324    575.29E0,
 325    576.00E0,
 326    625.55E0,
 327     20.15E0,
 328     28.78E0,
 329     29.57E0,
 330     37.41E0,
 331     39.12E0,
 332     50.24E0,
 333     61.38E0,
 334     66.25E0,
 335     73.42E0,
 336     95.52E0,
 337    107.32E0,
 338    122.04E0,
 339    134.03E0,
 340    163.19E0,
 341    163.48E0,
 342    175.70E0,
 343    179.86E0,
 344    211.27E0,
 345    217.78E0,
 346    219.14E0,
 347    262.52E0,
 348    268.01E0,
 349    268.62E0,
 350    336.25E0,
 351    337.23E0,
 352    339.33E0,
 353    427.38E0,
 354    428.58E0,
 355    432.68E0,
 356    528.99E0,
 357    531.08E0,
 358    628.34E0,
 359    253.24E0,
 360    273.13E0,
 361    273.66E0,
 362    282.10E0,
 363    346.62E0,
 364    347.19E0,
 365    348.78E0,
 366    351.18E0,
 367    450.10E0,
 368    450.35E0,
 369    451.92E0,
 370    455.56E0,
 371    552.22E0,
 372    553.56E0,
 373    555.74E0,
 374    652.59E0,
 375    656.20E0,
 376     14.13E0,
 377     20.41E0,
 378     31.30E0,
 379     33.84E0,
 380     39.70E0,
 381     48.83E0,
 382     54.50E0,
 383     60.41E0,
 384     72.77E0,
 385     75.25E0,
 386     86.84E0,
 387     94.88E0,
 388     96.40E0,
 389    117.37E0,
 390    139.08E0,
 391    147.73E0,
 392    158.63E0,
 393    161.84E0,
 394    192.11E0,
 395    206.76E0,
 396    209.07E0,
 397    213.32E0,
 398    226.44E0,
 399    237.12E0,
 400    330.90E0,
 401    358.72E0,
 402    370.77E0,
 403    372.72E0,
 404    396.24E0,
 405    416.59E0,
 406    484.02E0,
 407    495.47E0,
 408    514.78E0,
 409    515.65E0,
 410    519.47E0,
 411    544.47E0,
 412    560.11E0,
 413    620.77E0,
 414     18.97E0,
 415     28.93E0,
 416     33.91E0,
 417     40.03E0,
 418     44.66E0,
 419     49.87E0,
 420     55.16E0,
 421     60.90E0,
 422     72.08E0,
 423     85.15E0,
 424     97.06E0,
 425    119.63E0,
 426    133.27E0,
 427    143.84E0,
 428    161.91E0,
 429    180.67E0,
 430    198.44E0,
 431    226.86E0,
 432    229.65E0,
 433    258.27E0,
 434    273.77E0,
 435    339.15E0,
 436    350.13E0,
 437    362.75E0,
 438    371.03E0,
 439    393.32E0,
 440    448.53E0,
 441    473.78E0,
 442    511.12E0,
 443    524.70E0,
 444    548.75E0,
 445    551.64E0,
 446    574.02E0,
 447    623.86E0,
 448     21.46E0,
 449     24.33E0,
 450     33.43E0,
 451     39.22E0,
 452     44.18E0,
 453     55.02E0,
 454     94.33E0,
 455     96.44E0,
 456    118.82E0,
 457    128.48E0,
 458    141.94E0,
 459    156.92E0,
 460    171.65E0,
 461    190.00E0,
 462    223.26E0,
 463    223.88E0,
 464    231.50E0,
 465    265.05E0,
 466    269.44E0,
 467    271.78E0,
 468    273.46E0,
 469    334.61E0,
 470    339.79E0,
 471    349.52E0,
 472    358.18E0,
 473    377.98E0,
 474    394.77E0,
 475    429.66E0,
 476    468.22E0,
 477    487.27E0,
 478    519.54E0,
 479    523.03E0,
 480    612.99E0,
 481    638.59E0,
 482    641.36E0,
 483    622.05E0,
 484    631.50E0,
 485    663.97E0,
 486    646.9E0 ,
 487    748.29E0,
 488    749.21E0,
 489    750.14E0,
 490    647.04E0,
 491    646.89E0,
 492    746.9E0 ,
 493    748.43E0,
 494    747.35E0,
 495    749.27E0,
 496    647.61E0,
 497    747.78E0,
 498    750.51E0,
 499    851.37E0,
 500    845.97E0,
 501    847.54E0,
 502    849.93E0,
 503    851.61E0,
 504    849.75E0,
 505    850.98E0,
 506    848.23E0
 507 };
 508
 509
 510 int
 511 hahn1_f (const gsl_vector * x, void *params, gsl_vector * f)
 512 {
 513   double b[7];
 514   size_t i;
 515
 516   for (i = 0; i < 7; i++)
 517     {
 518       b[i] = gsl_vector_get(x, i);
 519     }
 520
 521   for (i = 0; i < 236; i++)
 522     {
 523       double x = hahn1_F0[i];
 524       double y = ((b[0] + x* (b[1]  + x * (b[2] + x * b[3])))
 525                   / (1 + x*(b[4]  + x *(b[5] + x*b[6]))));
 526       gsl_vector_set (f, i, hahn1_F1[i] - y);
 527     }
 528
 529   return GSL_SUCCESS;
 530 }
 531
 532 int
 533 hahn1_df (const gsl_vector * x, void *params, gsl_matrix * df)
 534 {
 535   double b[7];
 536   size_t i;
 537
 538   for (i = 0; i < 7; i++)
 539     {
 540       b[i] = gsl_vector_get(x, i);
 541     }
 542
 543   for (i = 0; i < 236; i++)
 544     {
 545       double x = hahn1_F0[i];
 546       double u = (b[0] + x*(b[1] + x*(b[2] + x * b[3])));
 547       double v = (1 + x*(b[4] + x*(b[5] + x*b[6])));
 548       gsl_matrix_set (df, i, 0, -1/v);
 549       gsl_matrix_set (df, i, 1, -x/v);
 550       gsl_matrix_set (df, i, 2, -x*x/v);
 551       gsl_matrix_set (df, i, 3, -x*x*x/v);
 552       gsl_matrix_set (df, i, 4, x*u/(v*v));
 553       gsl_matrix_set (df, i, 5, x*x*u/(v*v));
 554       gsl_matrix_set (df, i, 6, x*x*x*u/(v*v));
 555     }
 556
 557   return GSL_SUCCESS;
 558 }
 559
 560 int
 561 hahn1_fdf (const gsl_vector * x, void *params,
 562            gsl_vector * f, gsl_matrix * df)
 563 {
 564   hahn1_f (x, params, f);
 565   hahn1_df (x, params, df);
 566
 567   return GSL_SUCCESS;
 568 }