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

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