expressions/helpers.c

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