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

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