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

   1 const size_t kirby2_N = 151;
   2 const size_t kirby2_P = 5;
   3
   4 /* double kirby2_x0[5] = { 2, -0.1, 0.003, -0.001, 0.00001 }; */
   5
   6 double kirby2_x0[5] = { 1.5, -0.15, 0.0025, -0.0015, 0.00002 };
   7
   8 double kirby2_x[5] = {
   9   1.6745063063E+00,
  10   -1.3927397867E-01,
  11   2.5961181191E-03,
  12   -1.7241811870E-03,
  13   2.1664802578E-05
  14 };
  15
  16 double kirby2_sumsq = 3.9050739624E+00;
  17
  18 double kirby2_sigma[5] = {
  19   8.7989634338E-02,
  20   4.1182041386E-03,
  21   4.1856520458E-05,
  22   5.8931897355E-05,
  23   2.0129761919E-07
  24 };
  25
  26 double kirby2_F1[151] = {
  27        0.0082E0,
  28        0.0112E0,
  29        0.0149E0,
  30        0.0198E0,
  31        0.0248E0,
  32        0.0324E0,
  33        0.0420E0,
  34        0.0549E0,
  35        0.0719E0,
  36        0.0963E0,
  37        0.1291E0,
  38        0.1710E0,
  39        0.2314E0,
  40        0.3227E0,
  41        0.4809E0,
  42        0.7084E0,
  43        1.0220E0,
  44        1.4580E0,
  45        1.9520E0,
  46        2.5410E0,
  47        3.2230E0,
  48        3.9990E0,
  49        4.8520E0,
  50        5.7320E0,
  51        6.7270E0,
  52        7.8350E0,
  53        9.0250E0,
  54       10.2670E0,
  55       11.5780E0,
  56       12.9440E0,
  57       14.3770E0,
  58       15.8560E0,
  59       17.3310E0,
  60       18.8850E0,
  61       20.5750E0,
  62       22.3200E0,
  63       22.3030E0,
  64       23.4600E0,
  65       24.0600E0,
  66       25.2720E0,
  67       25.8530E0,
  68       27.1100E0,
  69       27.6580E0,
  70       28.9240E0,
  71       29.5110E0,
  72       30.7100E0,
  73       31.3500E0,
  74       32.5200E0,
  75       33.2300E0,
  76       34.3300E0,
  77       35.0600E0,
  78       36.1700E0,
  79       36.8400E0,
  80       38.0100E0,
  81       38.6700E0,
  82       39.8700E0,
  83       40.0300E0,
  84       40.5000E0,
  85       41.3700E0,
  86       41.6700E0,
  87       42.3100E0,
  88       42.7300E0,
  89       43.4600E0,
  90       44.1400E0,
  91       44.5500E0,
  92       45.2200E0,
  93       45.9200E0,
  94       46.3000E0,
  95       47.0000E0,
  96       47.6800E0,
  97       48.0600E0,
  98       48.7400E0,
  99       49.4100E0,
 100       49.7600E0,
 101       50.4300E0,
 102       51.1100E0,
 103       51.5000E0,
 104       52.1200E0,
 105       52.7600E0,
 106       53.1800E0,
 107       53.7800E0,
 108       54.4600E0,
 109       54.8300E0,
 110       55.4000E0,
 111       56.4300E0,
 112       57.0300E0,
 113       58.0000E0,
 114       58.6100E0,
 115       59.5800E0,
 116       60.1100E0,
 117       61.1000E0,
 118       61.6500E0,
 119       62.5900E0,
 120       63.1200E0,
 121       64.0300E0,
 122       64.6200E0,
 123       65.4900E0,
 124       66.0300E0,
 125       66.8900E0,
 126       67.4200E0,
 127       68.2300E0,
 128       68.7700E0,
 129       69.5900E0,
 130       70.1100E0,
 131       70.8600E0,
 132       71.4300E0,
 133       72.1600E0,
 134       72.7000E0,
 135       73.4000E0,
 136       73.9300E0,
 137       74.6000E0,
 138       75.1600E0,
 139       75.8200E0,
 140       76.3400E0,
 141       76.9800E0,
 142       77.4800E0,
 143       78.0800E0,
 144       78.6000E0,
 145       79.1700E0,
 146       79.6200E0,
 147       79.8800E0,
 148       80.1900E0,
 149       80.6600E0,
 150       81.2200E0,
 151       81.6600E0,
 152       82.1600E0,
 153       82.5900E0,
 154       83.1400E0,
 155       83.5000E0,
 156       84.0000E0,
 157       84.4000E0,
 158       84.8900E0,
 159       85.2600E0,
 160       85.7400E0,
 161       86.0700E0,
 162       86.5400E0,
 163       86.8900E0,
 164       87.3200E0,
 165       87.6500E0,
 166       88.1000E0,
 167       88.4300E0,
 168       88.8300E0,
 169       89.1200E0,
 170       89.5400E0,
 171       89.8500E0,
 172       90.2500E0,
 173       90.5500E0,
 174       90.9300E0,
 175       91.2000E0,
 176       91.5500E0,
 177       92.2000E0
 178 };
 179
 180
 181 double kirby2_F0[151] = {
 182     9.65E0,
 183    10.74E0,
 184    11.81E0,
 185    12.88E0,
 186    14.06E0,
 187    15.28E0,
 188    16.63E0,
 189    18.19E0,
 190    19.88E0,
 191    21.84E0,
 192    24.00E0,
 193    26.25E0,
 194    28.86E0,
 195    31.85E0,
 196    35.79E0,
 197    40.18E0,
 198    44.74E0,
 199    49.53E0,
 200    53.94E0,
 201    58.29E0,
 202    62.63E0,
 203    67.03E0,
 204    71.25E0,
 205    75.22E0,
 206    79.33E0,
 207    83.56E0,
 208    87.75E0,
 209    91.93E0,
 210    96.10E0,
 211   100.28E0,
 212   104.46E0,
 213   108.66E0,
 214   112.71E0,
 215   116.88E0,
 216   121.33E0,
 217   125.79E0,
 218   125.79E0,
 219   128.74E0,
 220   130.27E0,
 221   133.33E0,
 222   134.79E0,
 223   137.93E0,
 224   139.33E0,
 225   142.46E0,
 226   143.90E0,
 227   146.91E0,
 228   148.51E0,
 229   151.41E0,
 230   153.17E0,
 231   155.97E0,
 232   157.76E0,
 233   160.56E0,
 234   162.30E0,
 235   165.21E0,
 236   166.90E0,
 237   169.92E0,
 238   170.32E0,
 239   171.54E0,
 240   173.79E0,
 241   174.57E0,
 242   176.25E0,
 243   177.34E0,
 244   179.19E0,
 245   181.02E0,
 246   182.08E0,
 247   183.88E0,
 248   185.75E0,
 249   186.80E0,
 250   188.63E0,
 251   190.45E0,
 252   191.48E0,
 253   193.35E0,
 254   195.22E0,
 255   196.23E0,
 256   198.05E0,
 257   199.97E0,
 258   201.06E0,
 259   202.83E0,
 260   204.69E0,
 261   205.86E0,
 262   207.58E0,
 263   209.50E0,
 264   210.65E0,
 265   212.33E0,
 266   215.43E0,
 267   217.16E0,
 268   220.21E0,
 269   221.98E0,
 270   225.06E0,
 271   226.79E0,
 272   229.92E0,
 273   231.69E0,
 274   234.77E0,
 275   236.60E0,
 276   239.63E0,
 277   241.50E0,
 278   244.48E0,
 279   246.40E0,
 280   249.35E0,
 281   251.32E0,
 282   254.22E0,
 283   256.24E0,
 284   259.11E0,
 285   261.18E0,
 286   264.02E0,
 287   266.13E0,
 288   268.94E0,
 289   271.09E0,
 290   273.87E0,
 291   276.08E0,
 292   278.83E0,
 293   281.08E0,
 294   283.81E0,
 295   286.11E0,
 296   288.81E0,
 297   291.08E0,
 298   293.75E0,
 299   295.99E0,
 300   298.64E0,
 301   300.84E0,
 302   302.02E0,
 303   303.48E0,
 304   305.65E0,
 305   308.27E0,
 306   310.41E0,
 307   313.01E0,
 308   315.12E0,
 309   317.71E0,
 310   319.79E0,
 311   322.36E0,
 312   324.42E0,
 313   326.98E0,
 314   329.01E0,
 315   331.56E0,
 316   333.56E0,
 317   336.10E0,
 318   338.08E0,
 319   340.60E0,
 320   342.57E0,
 321   345.08E0,
 322   347.02E0,
 323   349.52E0,
 324   351.44E0,
 325   353.93E0,
 326   355.83E0,
 327   358.32E0,
 328   360.20E0,
 329   362.67E0,
 330   364.53E0,
 331   367.00E0,
 332   371.30E0
 333 };
 334
 335
 336 int
 337 kirby2_f (const gsl_vector * x, void *params, gsl_vector * f)
 338 {
 339   double b[5];
 340   size_t i;
 341
 342   for (i = 0; i < 5; i++)
 343     {
 344       b[i] = gsl_vector_get(x, i);
 345     }
 346
 347   for (i = 0; i < 151; i++)
 348     {
 349       double x = kirby2_F0[i];
 350       double y = ((b[0] + x* (b[1]  + x * b[2]))
 351                   / (1 + x*(b[3]  + x *b[4])));
 352       gsl_vector_set (f, i, kirby2_F1[i] - y);
 353     }
 354
 355   return GSL_SUCCESS;
 356 }
 357
 358 int
 359 kirby2_df (const gsl_vector * x, void *params, gsl_matrix * df)
 360 {
 361   double b[5];
 362   size_t i;
 363
 364   for (i = 0; i < 5; i++)
 365     {
 366       b[i] = gsl_vector_get(x, i);
 367     }
 368
 369   for (i = 0; i < 151; i++)
 370     {
 371       double x = kirby2_F0[i];
 372       double u = (b[0] + x*(b[1] + x*b[2]));
 373       double v = (1 + x*(b[3] + x*b[4]));
 374       gsl_matrix_set (df, i, 0, -1/v);
 375       gsl_matrix_set (df, i, 1, -x/v);
 376       gsl_matrix_set (df, i, 2, -x*x/v);
 377       gsl_matrix_set (df, i, 3, x*u/(v*v));
 378       gsl_matrix_set (df, i, 4, x*x*u/(v*v));
 379     }
 380
 381   return GSL_SUCCESS;
 382 }
 383
 384 int
 385 kirby2_fdf (const gsl_vector * x, void *params,
 386            gsl_vector * f, gsl_matrix * df)
 387 {
 388   kirby2_f (x, params, f);
 389   kirby2_df (x, params, df);
 390
 391   return GSL_SUCCESS;
 392 }