pintos-os.org Git - pspp/blob - src/language/expressions/helpers.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2008, 2010 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18 #include "helpers.h"
  19 #include <gsl/gsl_roots.h>
  20 #include <gsl/gsl_sf.h>
  21 #include <libpspp/assertion.h>
  22 #include <libpspp/pool.h>
  23 #include "private.h"
  24
  25 const struct substring empty_string = {NULL, 0};
  26
  27 double
  28 expr_ymd_to_ofs (double year, double month, double day)
  29 {
  30   int y = year;
  31   int m = month;
  32   int d = day;
  33   char *error;
  34   double ofs;
  35
  36   if (y != year || m != month || d != day)
  37     {
  38       msg (SE, _("One of the arguments to a DATE function is not an integer.  "
  39                  "The result will be system-missing."));
  40       return SYSMIS;
  41     }
  42
  43   ofs = calendar_gregorian_to_offset (y, m, d, &error);
  44   if (error != NULL)
  45     {
  46       msg (SE, "%s", error);
  47       free (error);
  48     }
  49   return ofs;
  50 }
  51
  52 double
  53 expr_ymd_to_date (double year, double month, double day)
  54 {
  55   double ofs = expr_ymd_to_ofs (year, month, day);
  56   return ofs != SYSMIS ? ofs * DAY_S : SYSMIS;
  57 }
  58
  59 double
  60 expr_wkyr_to_date (double week, double year)
  61 {
  62   int w = week;
  63
  64   if (w != week)
  65     {
  66       msg (SE, _("The week argument to DATE.WKYR is not an integer.  "
  67                  "The result will be system-missing."));
  68       return SYSMIS;
  69     }
  70   else if (w < 1 || w > 53)
  71     {
  72       msg (SE, _("The week argument to DATE.WKYR is outside the acceptable "
  73                  "range of 1 to 53.  "
  74                  "The result will be system-missing."));
  75       return SYSMIS;
  76     }
  77   else
  78     {
  79       double yr_1_1 = expr_ymd_to_ofs (year, 1, 1);
  80       if (yr_1_1 != SYSMIS)
  81         return DAY_S * (yr_1_1 + WEEK_DAY * (w - 1));
  82       else
  83         return SYSMIS;
  84     }
  85 }
  86
  87 double
  88 expr_yrday_to_date (double year, double yday)
  89 {
  90   int yd = yday;
  91
  92   if (yd != yday)
  93     {
  94       msg (SE, _("The day argument to DATE.YRDAY is not an integer.  "
  95                  "The result will be system-missing."));
  96       return SYSMIS;
  97     }
  98   else if (yd < 1 || yd > 366)
  99     {
 100       msg (SE, _("The day argument to DATE.YRDAY is outside the acceptable "
 101                  "range of 1 to 366.  "
 102                  "The result will be system-missing."));
 103       return SYSMIS;
 104     }
 105   else
 106     {
 107       double yr_1_1 = expr_ymd_to_ofs (year, 1, 1);
 108       if (yr_1_1 != SYSMIS)
 109         return DAY_S * (yr_1_1 + yd - 1.);
 110       else
 111         return SYSMIS;
 112     }
 113 }
 114
 115 double
 116 expr_yrmoda (double year, double month, double day)
 117 {
 118   if (year >= 0 && year <= 99)
 119     year += 1900;
 120   else if (year != (int) year && year > 47516)
 121     {
 122       msg (SE, _("The year argument to YRMODA is greater than 47516.  "
 123                  "The result will be system-missing."));
 124       return SYSMIS;
 125     }
 126
 127   return expr_ymd_to_ofs (year, month, day);
 128 }
 129 \f
 130 /* A date unit. */
 131 enum date_unit
 132   {
 133     DATE_YEARS,
 134     DATE_QUARTERS,
 135     DATE_MONTHS,
 136     DATE_WEEKS,
 137     DATE_DAYS,
 138     DATE_HOURS,
 139     DATE_MINUTES,
 140     DATE_SECONDS
 141   };
 142
 143 /* Stores in *UNIT the unit whose name is NAME.
 144    Return success. */
 145 static enum date_unit
 146 recognize_unit (struct substring name, enum date_unit *unit)
 147 {
 148   struct unit_name
 149     {
 150       enum date_unit unit;
 151       const struct substring name;
 152     };
 153   static const struct unit_name unit_names[] =
 154     {
 155       { DATE_YEARS, SS_LITERAL_INITIALIZER ("years") },
 156       { DATE_QUARTERS, SS_LITERAL_INITIALIZER ("quarters") },
 157       { DATE_MONTHS, SS_LITERAL_INITIALIZER ("months") },
 158       { DATE_WEEKS, SS_LITERAL_INITIALIZER ("weeks") },
 159       { DATE_DAYS, SS_LITERAL_INITIALIZER ("days") },
 160       { DATE_HOURS, SS_LITERAL_INITIALIZER ("hours") },
 161       { DATE_MINUTES, SS_LITERAL_INITIALIZER ("minutes") },
 162       { DATE_SECONDS, SS_LITERAL_INITIALIZER ("seconds") },
 163     };
 164   const int unit_name_cnt = sizeof unit_names / sizeof *unit_names;
 165
 166   const struct unit_name *un;
 167
 168   for (un = unit_names; un < &unit_names[unit_name_cnt]; un++)
 169     if (ss_equals_case (un->name, name))
 170       {
 171         *unit = un->unit;
 172         return true;
 173       }
 174
 175   /* TRANSLATORS: Don't translate the the actual unit names `weeks', `days' etc
 176         They must remain in their original English. */
 177   msg (SE, _("Unrecognized date unit `%.*s'.  "
 178              "Valid date units are `years', `quarters', `months', "
 179              "`weeks', `days', `hours', `minutes', and `seconds'."),
 180        (int) ss_length (name), ss_data (name));
 181   return false;
 182 }
 183
 184 /* Returns the number of whole years from DATE1 to DATE2,
 185    where a year is defined as the same or later month, day, and
 186    time of day. */
 187 static int
 188 year_diff (double date1, double date2)
 189 {
 190   int y1, m1, d1, yd1;
 191   int y2, m2, d2, yd2;
 192   int diff;
 193
 194   assert (date2 >= date1);
 195   calendar_offset_to_gregorian (date1 / DAY_S, &y1, &m1, &d1, &yd1);
 196   calendar_offset_to_gregorian (date2 / DAY_S, &y2, &m2, &d2, &yd2);
 197
 198   diff = y2 - y1;
 199   if (diff > 0)
 200     {
 201       int yd1 = 32 * m1 + d1;
 202       int yd2 = 32 * m2 + d2;
 203       if (yd2 < yd1
 204           || (yd2 == yd1 && fmod (date2, DAY_S) < fmod (date1, DAY_S)))
 205         diff--;
 206     }
 207   return diff;
 208 }
 209
 210 /* Returns the number of whole months from DATE1 to DATE2,
 211    where a month is defined as the same or later day and time of
 212    day. */
 213 static int
 214 month_diff (double date1, double date2)
 215 {
 216   int y1, m1, d1, yd1;
 217   int y2, m2, d2, yd2;
 218   int diff;
 219
 220   assert (date2 >= date1);
 221   calendar_offset_to_gregorian (date1 / DAY_S, &y1, &m1, &d1, &yd1);
 222   calendar_offset_to_gregorian (date2 / DAY_S, &y2, &m2, &d2, &yd2);
 223
 224   diff = ((y2 * 12) + m2) - ((y1 * 12) + m1);
 225   if (diff > 0
 226       && (d2 < d1
 227           || (d2 == d1 && fmod (date2, DAY_S) < fmod (date1, DAY_S))))
 228     diff--;
 229   return diff;
 230 }
 231
 232 /* Returns the number of whole quarter from DATE1 to DATE2,
 233    where a quarter is defined as three months. */
 234 static int
 235 quarter_diff (double date1, double date2)
 236 {
 237   return month_diff (date1, date2) / 3;
 238 }
 239
 240 /* Returns the number of seconds in the given UNIT. */
 241 static int
 242 date_unit_duration (enum date_unit unit)
 243 {
 244   switch (unit)
 245     {
 246     case DATE_WEEKS:
 247       return WEEK_S;
 248
 249     case DATE_DAYS:
 250       return DAY_S;
 251
 252     case DATE_HOURS:
 253       return H_S;
 254
 255     case DATE_MINUTES:
 256       return MIN_S;
 257
 258     case DATE_SECONDS:
 259       return 1;
 260
 261     default:
 262       NOT_REACHED ();
 263     }
 264 }
 265
 266 /* Returns the span from DATE1 to DATE2 in terms of UNIT_NAME. */
 267 double
 268 expr_date_difference (double date1, double date2, struct substring unit_name)
 269 {
 270   enum date_unit unit;
 271
 272   if (!recognize_unit (unit_name, &unit))
 273     return SYSMIS;
 274
 275   switch (unit)
 276     {
 277     case DATE_YEARS:
 278       return (date2 >= date1
 279               ? year_diff (date1, date2)
 280               : -year_diff (date2, date1));
 281
 282     case DATE_QUARTERS:
 283       return (date2 >= date1
 284               ? quarter_diff (date1, date2)
 285               : -quarter_diff (date2, date1));
 286
 287     case DATE_MONTHS:
 288       return (date2 >= date1
 289               ? month_diff (date1, date2)
 290               : -month_diff (date2, date1));
 291
 292     case DATE_WEEKS:
 293     case DATE_DAYS:
 294     case DATE_HOURS:
 295     case DATE_MINUTES:
 296     case DATE_SECONDS:
 297       return trunc ((date2 - date1) / date_unit_duration (unit));
 298     }
 299
 300   NOT_REACHED ();
 301 }
 302
 303 /* How to deal with days out of range for a given month. */
 304 enum date_sum_method
 305   {
 306     SUM_ROLLOVER,       /* Roll them over to the next month. */
 307     SUM_CLOSEST         /* Use the last day of the month. */
 308   };
 309
 310 /* Stores in *METHOD the method whose name is NAME.
 311    Return success. */
 312 static bool
 313 recognize_method (struct substring method_name, enum date_sum_method *method)
 314 {
 315   if (ss_equals_case (method_name, ss_cstr ("closest")))
 316     {
 317       *method = SUM_CLOSEST;
 318       return true;
 319     }
 320   else if (ss_equals_case (method_name, ss_cstr ("rollover")))
 321     {
 322       *method = SUM_ROLLOVER;
 323       return true;
 324     }
 325   else
 326     {
 327       msg (SE, _("Invalid DATESUM method.  "
 328                  "Valid choices are `closest' and `rollover'."));
 329       return false;
 330     }
 331 }
 332
 333 /* Returns DATE advanced by the given number of MONTHS, with
 334    day-of-month overflow resolved using METHOD. */
 335 static double
 336 add_months (double date, int months, enum date_sum_method method)
 337 {
 338   int y, m, d, yd;
 339   double output;
 340   char *error;
 341
 342   calendar_offset_to_gregorian (date / DAY_S, &y, &m, &d, &yd);
 343   y += months / 12;
 344   m += months % 12;
 345   if (m < 1)
 346     {
 347       m += 12;
 348       y--;
 349     }
 350   else if (m > 12)
 351     {
 352       m -= 12;
 353       y++;
 354     }
 355   assert (m >= 1 && m <= 12);
 356
 357   if (method == SUM_CLOSEST && d > calendar_days_in_month (y, m))
 358     d = calendar_days_in_month (y, m);
 359
 360   output = calendar_gregorian_to_offset (y, m, d, &error);
 361   if (output != SYSMIS)
 362     output = (output * DAY_S) + fmod (date, DAY_S);
 363   else
 364     {
 365       msg (SE, "%s", error);
 366       free (error);
 367     }
 368   return output;
 369 }
 370
 371 /* Returns DATE advanced by the given QUANTITY of units given in
 372    UNIT_NAME, with day-of-month overflow resolved using
 373    METHOD_NAME. */
 374 double
 375 expr_date_sum (double date, double quantity, struct substring unit_name,
 376                struct substring method_name)
 377 {
 378   enum date_unit unit;
 379   enum date_sum_method method;
 380
 381   if (!recognize_unit (unit_name, &unit)
 382       || !recognize_method (method_name, &method))
 383     return SYSMIS;
 384
 385   switch (unit)
 386     {
 387     case DATE_YEARS:
 388       return add_months (date, trunc (quantity) * 12, method);
 389
 390     case DATE_QUARTERS:
 391       return add_months (date, trunc (quantity) * 3, method);
 392
 393     case DATE_MONTHS:
 394       return add_months (date, trunc (quantity), method);
 395
 396     case DATE_WEEKS:
 397     case DATE_DAYS:
 398     case DATE_HOURS:
 399     case DATE_MINUTES:
 400     case DATE_SECONDS:
 401       return date + quantity * date_unit_duration (unit);
 402     }
 403
 404   NOT_REACHED ();
 405 }
 406
 407 int
 408 compare_string_3way (const struct substring *a, const struct substring *b)
 409 {
 410   size_t i;
 411
 412   for (i = 0; i < a->length && i < b->length; i++)
 413     if (a->string[i] != b->string[i])
 414       return a->string[i] < b->string[i] ? -1 : 1;
 415   for (; i < a->length; i++)
 416     if (a->string[i] != ' ')
 417       return 1;
 418   for (; i < b->length; i++)
 419     if (b->string[i] != ' ')
 420       return -1;
 421   return 0;
 422 }
 423
 424 size_t
 425 count_valid (double *d, size_t d_cnt)
 426 {
 427   size_t valid_cnt;
 428   size_t i;
 429
 430   valid_cnt = 0;
 431   for (i = 0; i < d_cnt; i++)
 432     valid_cnt += is_valid (d[i]);
 433   return valid_cnt;
 434 }
 435
 436 struct substring
 437 alloc_string (struct expression *e, size_t length)
 438 {
 439   struct substring s;
 440   s.length = length;
 441   s.string = pool_alloc (e->eval_pool, length);
 442   return s;
 443 }
 444
 445 struct substring
 446 copy_string (struct expression *e, const char *old, size_t length)
 447 {
 448   struct substring s = alloc_string (e, length);
 449   memcpy (s.string, old, length);
 450   return s;
 451 }
 452
 453 /* Returns the noncentral beta cumulative distribution function
 454    value for the given arguments.
 455
 456    FIXME: The accuracy of this function is not entirely
 457    satisfactory.  We only match the example values given in AS
 458    310 to the first 5 significant digits. */
 459 double
 460 ncdf_beta (double x, double a, double b, double lambda)
 461 {
 462   double c;
 463
 464   if (x <= 0. || x >= 1. || a <= 0. || b <= 0. || lambda <= 0.)
 465     return SYSMIS;
 466
 467   c = lambda / 2.;
 468   if (lambda < 54.)
 469     {
 470       /* Algorithm AS 226. */
 471       double x0, a0, beta, temp, gx, q, ax, sumq, sum;
 472       double err_max = 2 * DBL_EPSILON;
 473       double err_bound;
 474       int iter_max = 100;
 475       int iter;
 476
 477       x0 = floor (c - 5.0 * sqrt (c));
 478       if (x0 < 0.)
 479         x0 = 0.;
 480       a0 = a + x0;
 481       beta = (gsl_sf_lngamma (a0)
 482               + gsl_sf_lngamma (b)
 483               - gsl_sf_lngamma (a0 + b));
 484       temp = gsl_sf_beta_inc (a0, b, x);
 485       gx = exp (a0 * log (x) + b * log (1. - x) - beta - log (a0));
 486       if (a0 >= a)
 487         q = exp (-c + x0 * log (c)) - gsl_sf_lngamma (x0 + 1.);
 488       else
 489         q = exp (-c);
 490       ax = q * temp;
 491       sumq = 1. - q;
 492       sum = ax;
 493
 494       iter = 0;
 495       do
 496         {
 497           iter++;
 498           temp -= gx;
 499           gx = x * (a + b + iter - 1.) * gx / (a + iter);
 500           q *= c / iter;
 501           sumq -= q;
 502           ax = temp * q;
 503           sum += ax;
 504
 505           err_bound = (temp - gx) * sumq;
 506         }
 507       while (iter < iter_max && err_bound > err_max);
 508
 509       return sum;
 510     }
 511   else
 512     {
 513       /* Algorithm AS 310. */
 514       double m, m_sqrt;
 515       int iter, iter_lower, iter_upper, iter1, iter2, j;
 516       double t, q, r, psum, beta, s1, gx, fx, temp, ftemp, t0, s0, sum, s;
 517       double err_bound;
 518       double err_max = 2 * DBL_EPSILON;
 519
 520       iter = 0;
 521
 522       m = floor (c + .5);
 523       m_sqrt = sqrt (m);
 524       iter_lower = m - 5. * m_sqrt;
 525       iter_upper = m + 5. * m_sqrt;
 526
 527       t = -c + m * log (c) - gsl_sf_lngamma (m + 1.);
 528       q = exp (t);
 529       r = q;
 530       psum = q;
 531       beta = (gsl_sf_lngamma (a + m)
 532               + gsl_sf_lngamma (b)
 533               - gsl_sf_lngamma (a + m + b));
 534       s1 = (a + m) * log (x) + b * log (1. - x) - log (a + m) - beta;
 535       fx = gx = exp (s1);
 536       ftemp = temp = gsl_sf_beta_inc (a + m, b, x);
 537       iter++;
 538       sum = q * temp;
 539       iter1 = m;
 540
 541       while (iter1 >= iter_lower && q >= err_max)
 542         {
 543           q = q * iter1 / c;
 544           iter++;
 545           gx = (a + iter1) / (x * (a + b + iter1 - 1.)) * gx;
 546           iter1--;
 547           temp += gx;
 548           psum += q;
 549           sum += q * temp;
 550         }
 551
 552       t0 = (gsl_sf_lngamma (a + b)
 553             - gsl_sf_lngamma (a + 1.)
 554             - gsl_sf_lngamma (b));
 555       s0 = a * log (x) + b * log (1. - x);
 556
 557       s = 0.;
 558       for (j = 0; j < iter1; j++)
 559         {
 560           double t1;
 561           s += exp (t0 + s0 + j * log (x));
 562           t1 = log (a + b + j) - log (a + 1. + j) + t0;
 563           t0 = t1;
 564         }
 565
 566       err_bound = (1. - gsl_sf_gamma_inc_P (iter1, c)) * (temp + s);
 567       q = r;
 568       temp = ftemp;
 569       gx = fx;
 570       iter2 = m;
 571       for (;;)
 572         {
 573           double ebd = err_bound + (1. - psum) * temp;
 574           if (ebd < err_max || iter >= iter_upper)
 575             break;
 576
 577           iter2++;
 578           iter++;
 579           q = q * c / iter2;
 580           psum += q;
 581           temp -= gx;
 582           gx = x * (a + b + iter2 - 1.) / (a + iter2) * gx;
 583           sum += q * temp;
 584         }
 585
 586       return sum;
 587     }
 588 }
 589
 590 double
 591 cdf_bvnor (double x0, double x1, double r)
 592 {
 593   double z = pow2 (x0) - 2. * r * x0 * x1 + pow2 (x1);
 594   return exp (-z / (2. * (1 - r * r))) * (2. * M_PI * sqrt (1 - r * r));
 595 }
 596
 597 double
 598 idf_fdist (double P, double df1, double df2)
 599 {
 600   double temp = gsl_cdf_beta_Pinv (P, df1 / 2, df2 / 2);
 601   return temp * df2 / ((1. - temp) * df1);
 602 }
 603
 604 /*
 605  *  Mathlib : A C Library of Special Functions
 606  *  Copyright (C) 1998 Ross Ihaka
 607  *  Copyright (C) 2000 The R Development Core Team
 608  *
 609  *  This program is free software; you can redistribute it and/or
 610  *  modify
 611  *  it under the terms of the GNU General Public License as
 612  *  published by
 613  *  the Free Software Foundation; either version 2 of the
 614  *  License, or
 615  *  (at your option) any later version.
 616  *
 617  *  This program is distributed in the hope that it will be
 618  *  useful,
 619  *  but WITHOUT ANY WARRANTY; without even the implied warranty
 620  *  of
 621  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 622  *  GNU General Public License for more details.
 623  *
 624  *  You should have received a copy of the GNU General Public
 625  *  License
 626  *  along with this program; if not, write to the Free Software
 627  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 628  *  02110-1301 USA.
 629  */
 630
 631 /* Returns the density of the noncentral beta distribution with
 632    noncentrality parameter LAMBDA. */
 633 double
 634 npdf_beta (double x, double a, double b, double lambda)
 635 {
 636   if (lambda < 0. || a <= 0. || b <= 0.)
 637     return SYSMIS;
 638   else if (lambda == 0.)
 639     return gsl_ran_beta_pdf (x, a, b);
 640   else
 641     {
 642       double max_error = 2 * DBL_EPSILON;
 643       int max_iter = 200;
 644       double term = gsl_ran_beta_pdf (x, a, b);
 645       double lambda2 = 0.5 * lambda;
 646       double weight = exp (-lambda2);
 647       double sum = weight * term;
 648       double psum = weight;
 649       int k;
 650       for (k = 1; k <= max_iter && 1 - psum < max_error; k++) {
 651         weight *= lambda2 / k;
 652         term *= x * (a + b) / a;
 653         sum += weight * term;
 654         psum += weight;
 655         a += 1;
 656       }
 657       return sum;
 658     }
 659 }