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

   1 const size_t nelson_N = 128;
   2 const size_t nelson_P = 3;
   3
   4 /* double nelson_x0[3] = { 2, 0.0001, -0.01}; */
   5
   6 double nelson_x0[3] = { 2.5 , 0.000000005, -0.01};
   7
   8 double nelson_x[3] = {
   9   2.5906836021E+00,
  10   5.6177717026E-09,
  11   -5.7701013174E-02
  12 };
  13
  14 double nelson_sumsq = 3.7976833176E+00;
  15
  16 double nelson_sigma[3] = {
  17   1.9149996413E-02,
  18   6.1124096540E-09,
  19   3.9572366543E-03
  20 };
  21
  22 double nelson_F[128][3] = {
  23   {  15.00E0,      1E0,      180E0},
  24   {  17.00E0,      1E0,      180E0},
  25   {  15.50E0,      1E0,      180E0},
  26   {  16.50E0,      1E0,      180E0},
  27   {  15.50E0,      1E0,      225E0},
  28   {  15.00E0,      1E0,      225E0},
  29   {  16.00E0,      1E0,      225E0},
  30   {  14.50E0,      1E0,      225E0},
  31   {  15.00E0,      1E0,      250E0},
  32   {  14.50E0,      1E0,      250E0},
  33   {  12.50E0,      1E0,      250E0},
  34   {  11.00E0,      1E0,      250E0},
  35   {  14.00E0,      1E0,      275E0},
  36   {  13.00E0,      1E0,      275E0},
  37   {  14.00E0,      1E0,      275E0},
  38   {  11.50E0,      1E0,      275E0},
  39   {  14.00E0,      2E0,      180E0},
  40   {  16.00E0,      2E0,      180E0},
  41   {  13.00E0,      2E0,      180E0},
  42   {  13.50E0,      2E0,      180E0},
  43   {  13.00E0,      2E0,      225E0},
  44   {  13.50E0,      2E0,      225E0},
  45   {  12.50E0,      2E0,      225E0},
  46   {  12.50E0,      2E0,      225E0},
  47   {  12.50E0,      2E0,      250E0},
  48   {  12.00E0,      2E0,      250E0},
  49   {  11.50E0,      2E0,      250E0},
  50   {  12.00E0,      2E0,      250E0},
  51   {  13.00E0,      2E0,      275E0},
  52   {  11.50E0,      2E0,      275E0},
  53   {  13.00E0,      2E0,      275E0},
  54   {  12.50E0,      2E0,      275E0},
  55   {  13.50E0,      4E0,      180E0},
  56   {  17.50E0,      4E0,      180E0},
  57   {  17.50E0,      4E0,      180E0},
  58   {  13.50E0,      4E0,      180E0},
  59   {  12.50E0,      4E0,      225E0},
  60   {  12.50E0,      4E0,      225E0},
  61   {  15.00E0,      4E0,      225E0},
  62   {  13.00E0,      4E0,      225E0},
  63   {  12.00E0,      4E0,      250E0},
  64   {  13.00E0,      4E0,      250E0},
  65   {  12.00E0,      4E0,      250E0},
  66   {  13.50E0,      4E0,      250E0},
  67   {  10.00E0,      4E0,      275E0},
  68   {  11.50E0,      4E0,      275E0},
  69   {  11.00E0,      4E0,      275E0},
  70   {   9.50E0,      4E0,      275E0},
  71   {  15.00E0,      8E0,      180E0},
  72   {  15.00E0,      8E0,      180E0},
  73   {  15.50E0,      8E0,      180E0},
  74   {  16.00E0,      8E0,      180E0},
  75   {  13.00E0,      8E0,      225E0},
  76   {  10.50E0,      8E0,      225E0},
  77   {  13.50E0,      8E0,      225E0},
  78   {  14.00E0,      8E0,      225E0},
  79   {  12.50E0,      8E0,      250E0},
  80   {  12.00E0,      8E0,      250E0},
  81   {  11.50E0,      8E0,      250E0},
  82   {  11.50E0,      8E0,      250E0},
  83   {   6.50E0,      8E0,      275E0},
  84   {   5.50E0,      8E0,      275E0},
  85   {   6.00E0,      8E0,      275E0},
  86   {   6.00E0,      8E0,      275E0},
  87   {  18.50E0,     16E0,      180E0},
  88   {  17.00E0,     16E0,      180E0},
  89   {  15.30E0,     16E0,      180E0},
  90   {  16.00E0,     16E0,      180E0},
  91   {  13.00E0,     16E0,      225E0},
  92   {  14.00E0,     16E0,      225E0},
  93   {  12.50E0,     16E0,      225E0},
  94   {  11.00E0,     16E0,      225E0},
  95   {  12.00E0,     16E0,      250E0},
  96   {  12.00E0,     16E0,      250E0},
  97   {  11.50E0,     16E0,      250E0},
  98   {  12.00E0,     16E0,      250E0},
  99   {   6.00E0,     16E0,      275E0},
 100   {   6.00E0,     16E0,      275E0},
 101   {   5.00E0,     16E0,      275E0},
 102   {   5.50E0,     16E0,      275E0},
 103   {  12.50E0,     32E0,      180E0},
 104   {  13.00E0,     32E0,      180E0},
 105   {  16.00E0,     32E0,      180E0},
 106   {  12.00E0,     32E0,      180E0},
 107   {  11.00E0,     32E0,      225E0},
 108   {   9.50E0,     32E0,      225E0},
 109   {  11.00E0,     32E0,      225E0},
 110   {  11.00E0,     32E0,      225E0},
 111   {  11.00E0,     32E0,      250E0},
 112   {  10.00E0,     32E0,      250E0},
 113   {  10.50E0,     32E0,      250E0},
 114   {  10.50E0,     32E0,      250E0},
 115   {   2.70E0,     32E0,      275E0},
 116   {   2.70E0,     32E0,      275E0},
 117   {   2.50E0,     32E0,      275E0},
 118   {   2.40E0,     32E0,      275E0},
 119   {  13.00E0,     48E0,      180E0},
 120   {  13.50E0,     48E0,      180E0},
 121   {  16.50E0,     48E0,      180E0},
 122   {  13.60E0,     48E0,      180E0},
 123   {  11.50E0,     48E0,      225E0},
 124   {  10.50E0,     48E0,      225E0},
 125   {  13.50E0,     48E0,      225E0},
 126   {  12.00E0,     48E0,      225E0},
 127   {   7.00E0,     48E0,      250E0},
 128   {   6.90E0,     48E0,      250E0},
 129   {   8.80E0,     48E0,      250E0},
 130   {   7.90E0,     48E0,      250E0},
 131   {   1.20E0,     48E0,      275E0},
 132   {   1.50E0,     48E0,      275E0},
 133   {   1.00E0,     48E0,      275E0},
 134   {   1.50E0,     48E0,      275E0},
 135   {  13.00E0,     64E0,      180E0},
 136   {  12.50E0,     64E0,      180E0},
 137   {  16.50E0,     64E0,      180E0},
 138   {  16.00E0,     64E0,      180E0},
 139   {  11.00E0,     64E0,      225E0},
 140   {  11.50E0,     64E0,      225E0},
 141   {  10.50E0,     64E0,      225E0},
 142   {  10.00E0,     64E0,      225E0},
 143   {   7.27E0,     64E0,      250E0},
 144   {   7.50E0,     64E0,      250E0},
 145   {   6.70E0,     64E0,      250E0},
 146   {   7.60E0,     64E0,      250E0},
 147   {   1.50E0,     64E0,      275E0},
 148   {   1.00E0,     64E0,      275E0},
 149   {   1.20E0,     64E0,      275E0},
 150   {   1.20E0,     64E0,      275E0}
 151 };
 152
 153
 154 int
 155 nelson_f (const gsl_vector * x, void *params, gsl_vector * f)
 156 {
 157   double b[3];
 158   size_t i;
 159
 160   for (i = 0; i < 3; i++)
 161     {
 162       b[i] = gsl_vector_get(x, i);
 163     }
 164
 165
 166
 167   for (i = 0; i < 128; i++)
 168     {
 169       double x1 = nelson_F[i][1];
 170       double x2 = nelson_F[i][2];
 171
 172       double y = b[0] - b[1] * x1 * (b[2]*x2 < -100) ? 0.0 : exp(-b[2] * x2);
 173
 174       gsl_vector_set (f, i, log(nelson_F[i][0]) - y);
 175     }
 176
 177   return GSL_SUCCESS;
 178 }
 179
 180 int
 181 nelson_df (const gsl_vector * x, void *params, gsl_matrix * df)
 182 {
 183   double b[3];
 184   size_t i;
 185
 186   for (i = 0; i < 3; i++)
 187     {
 188       b[i] = gsl_vector_get(x, i);
 189     }
 190
 191   for (i = 0; i < 128; i++)
 192     {
 193       double x1 = nelson_F[i][1];
 194       double x2 = nelson_F[i][2];
 195
 196       gsl_matrix_set (df, i, 0, -1.0);
 197       gsl_matrix_set (df, i, 1, x1 * exp(-b[2] * x2));
 198       gsl_matrix_set (df, i, 2, -b[1] * x1 * x2 * exp(-b[2] * x2));
 199     }
 200
 201   return GSL_SUCCESS;
 202 }
 203
 204 int
 205 nelson_fdf (const gsl_vector * x, void *params,
 206            gsl_vector * f, gsl_matrix * df)
 207 {
 208   nelson_f (x, params, f);
 209   nelson_df (x, params, df);
 210
 211   return GSL_SUCCESS;
 212 }
 213
 214
 215
 216
 217