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