src/nws_portability/Forecast/mse_forc.c

   1 #include <unistd.h>
   2
   3 #include <stdlib.h>
   4 #include <stdio.h>
   5 #include <math.h>
   6 #include <string.h>
   7 #include <strings.h>
   8
   9
  10
  11 #include "forc.h"
  12
  13 struct mse_state
  14 {
  15         forcl forc_list;
  16         int win;
  17         int method;
  18         int last_index;         /* optimization */
  19         double last_forecast;   /* optimization */
  20 };
  21
  22
  23 /*
  24  * forward refs for derived forecaster types
  25  */
  26 int TotalMSEForecast(char *state, double *forecast);
  27 int LocalMAEForecast(char *state, double *forecast);
  28
  29 int DerivedMethod(char *state);
  30
  31 /*
  32  * the MSE forecast chooses the best forecast from the derived
  33  * forecasters, one of which is the total MSE forecast
  34  */
  35 double
  36 MSEForecast(char *i_forcl)
  37 {
  38         int i;
  39         int min_i;
  40         double min_se;
  41         double sq_err;
  42         double forecast;
  43         forcl l_forcl;
  44         int ferr;
  45         fbuff series;
  46
  47         l_forcl = (forcl)i_forcl;
  48
  49         /*
  50          * first, check to see if there is any data -- the forecaster
  51          * routines may be called with an empty series
  52          */
  53         if(F_COUNT(l_forcl->forcs[0]->series) == 0)
  54         {
  55                 forecast = FORE_ERROR_VAL;
  56                 return(forecast);
  57         }
  58
  59
  60
  61         min_se = DBIG_VAL;
  62         min_i = 0;
  63
  64         /*
  65          * find the forecaster that has the lowest current
  66          * mean square error
  67          */
  68         for(i=0; i < l_forcl->derived_count; i++)
  69         {
  70                 series = l_forcl->derived_forcs[i]->se_series;
  71                 if(F_COUNT(series) > 0)
  72                 {
  73                         sq_err =
  74                         F_VAL(series,F_FIRST(series)) / F_COUNT(series);
  75                 }
  76                 else
  77                 {
  78                         sq_err = DBIG_VAL;
  79                 }
  80
  81                 if(sq_err < min_se)
  82                 {
  83                         min_i = i;
  84                         min_se = sq_err;
  85                 }
  86         }
  87
  88         /*
  89          * report the best
  90          */
  91         ferr =
  92 (l_forcl->derived_forcs[min_i])->forecast((l_forcl->derived_forcs[min_i])->state,
  93                                   &forecast);
  94
  95         /*
  96          * if an error has occurred, print and error and report
  97          * FORE_ERROR_VAL
  98          */
  99         if(ferr == 0)
 100         {
 101                 /*
 102                  * since we know there is at least one data item,
 103                  * return the last on error
 104                  */
 105                 forecast = F_VAL(l_forcl->forcs[0]->series,
 106                                 F_FIRST(l_forcl->forcs[0]->series));
 107                 return(forecast);
 108         }
 109
 110         /*
 111          * update the method that was used
 112          */
 113         l_forcl->total_mse_method =
 114                 DerivedMethod(l_forcl->derived_forcs[min_i]->state);
 115
 116         return(forecast);
 117 }
 118
 119 double
 120 MSEError(char *state)
 121 {
 122         forcl s;
 123         double val;
 124
 125         s = (forcl)state;
 126
 127         if(s->total_mse->count > 0.0)
 128         {
 129                 val = s->total_mse->se / s->total_mse->count;
 130         }
 131         else
 132         {
 133                 val = 0.0;
 134         }
 135
 136         return(val);
 137 }
 138
 139 int
 140 MSEMethod(char *state)
 141 {
 142         forcl s;
 143
 144         s = (forcl)state;
 145
 146         return(s->total_mse_method);
 147 }
 148
 149 /*
 150  * for forecaster interface
 151  */
 152 int
 153 TotalMSEForecast(char *state, double *forecast)
 154 {
 155         *forecast = MSEForecast(state);
 156
 157         if(*forecast == FORE_ERROR_VAL)
 158         {
 159                 return(0);
 160         }
 161         else
 162         {
 163                 return(1);
 164         }
 165 }
 166
 167 double
 168 MAEForecast(char *i_forcl)
 169 {
 170         int i;
 171         int min_i;
 172         double min_ae;
 173         double err;
 174         double forecast;
 175         forcl l_forcl;
 176         int ferr;
 177         fbuff series;
 178
 179         l_forcl = (forcl)i_forcl;
 180
 181         min_ae = DBIG_VAL;
 182         min_i = 0;
 183
 184         /*
 185          * find the forecaster that has the lowest current
 186          * mean absolute error
 187          */
 188         for(i=0; i < l_forcl->derived_count; i++)
 189         {
 190                 series = l_forcl->derived_forcs[i]->ae_series;
 191                 if(F_COUNT(series) > 0)
 192                 {
 193                         err =
 194                         F_VAL(series,F_FIRST(series)) / F_COUNT(series);
 195                 }
 196                 else
 197                 {
 198                         err = DBIG_VAL;
 199                 }
 200
 201                 if(err < min_ae)
 202                 {
 203                         min_i = i;
 204                         min_ae = err;
 205                 }
 206         }
 207
 208         /*
 209          * report the best
 210          */
 211         ferr =
 212 (l_forcl->derived_forcs[min_i])->forecast((l_forcl->derived_forcs[min_i])->state,
 213                                   &forecast);
 214
 215         /*
 216          * if an error has occurred, print and error and report
 217          * FORE_ERROR_VAL
 218          */
 219         if(ferr == 0)
 220         {
 221                 forecast = F_VAL(l_forcl->forcs[0]->series,
 222                                 F_FIRST(l_forcl->forcs[0]->series));
 223         }
 224
 225         /*
 226          * update the method that was used
 227          */
 228         l_forcl->total_mae_method =
 229                 DerivedMethod(l_forcl->derived_forcs[min_i]->state);
 230         return(forecast);
 231 }
 232
 233 double
 234 MAEError(char *state)
 235 {
 236         forcl s;
 237         double val;
 238
 239         s = (forcl)state;
 240
 241         if(s->total_mae->count > 0.0)
 242         {
 243                 val = s->total_mae->ae / s->total_mae->count;
 244         }
 245         else
 246         {
 247                 val = 0.0;
 248         }
 249
 250         return(val);
 251 }
 252
 253 int
 254 MAEMethod(char *state)
 255 {
 256         forcl s;
 257
 258         s = (forcl)state;
 259
 260         return(s->total_mae_method);
 261 }
 262
 263 int
 264 TotalMAEForecast(char *state, double *forecast)
 265 {
 266         *forecast = MAEForecast(state);
 267
 268         if(*forecast == FORE_ERROR_VAL)
 269         {
 270                 return(0);
 271         }
 272         else
 273         {
 274                 return(1);
 275         }
 276 }
 277
 278 char *
 279 InitWinMSE(fbuff series, fbuff time_stamps, char *params)
 280 {
 281         struct mse_state *state;
 282         int pcount;
 283         forcl l_forcl;
 284         int win;
 285
 286         state = (struct mse_state *)malloc(sizeof(struct mse_state));
 287
 288         if(state == NULL)
 289         {
 290                 return(NULL);
 291         }
 292
 293         /*
 294          * first parameter is the address of the forcl and the
 295          * second is the window size
 296          */
 297         pcount = sscanf(params,"%p %d",
 298                         &l_forcl,
 299                         &win);
 300
 301         if(pcount != 2)
 302         {
 303                 free(state);
 304                 return(NULL);
 305         }
 306
 307         state->forc_list = l_forcl;
 308         state->win = win;
 309
 310         return((char *)state);
 311 }
 312
 313 void
 314 FreeWinMSE(char *state)
 315 {
 316         free(state);
 317 }
 318
 319 char *
 320 InitWinMAE(fbuff series, fbuff time_stamps, char *params)
 321 {
 322         return(InitWinMSE(series,time_stamps,params));
 323 }
 324
 325 void
 326 FreeWinMAE(char *state)
 327 {
 328         FreeWinMSE(state);
 329 }
 330
 331 int
 332 LocalWinMSEForecast(char *state, double *out_f)
 333 {
 334         struct mse_state *s;
 335         int i;
 336         int win;
 337         int min_i;
 338         double min_se;
 339         double sq_err;
 340         double forecast;
 341         forcl l_forcl;
 342         int ferr;
 343         double temp_err;
 344         int windex;
 345         fbuff series;
 346
 347         s = (struct mse_state *)state;
 348
 349         l_forcl = s->forc_list;
 350         win = s->win;
 351
 352
 353         min_se = DBIG_VAL;
 354         min_i = 0;
 355
 356
 357         /*
 358          * find the forecaster that has the lowest current
 359          * mean square error --
 360          *
 361          * if the win size > 0, look only over the window
 362          * using the cumulative series.  Otherwise, use the total
 363          * cumulative value
 364          */
 365         for(i=0; i < l_forcl->count; i++)
 366         {
 367                 if((l_forcl->forcs[i])->count > 0.0)
 368                 {
 369                         series =
 370                         (l_forcl->forcs[i])->se_series;
 371
 372                         /*
 373                          * make sure we only window as much
 374                          * as there is
 375                          */
 376                         if(win > (F_COUNT(series) - 1))
 377                                 win = F_COUNT(series) - 1;
 378
 379                         if(win == 0)
 380                         {
 381                                 /*
 382                                  * get the first value which is
 383                                  * the current cumulative total
 384                                  */
 385                                 temp_err =
 386                                 F_VAL(series, F_FIRST(series));
 387                                 sq_err = temp_err /
 388                                         (l_forcl->forcs[i])->count;
 389                         }
 390                         else
 391                         {
 392                                 /*
 393                                  * get the index of the end of the
 394                                  * window
 395                                  */
 396                                 windex = MODPLUS(F_FIRST(series),
 397                                                  win,
 398                                                  F_SIZE(series));
 399                                 /*
 400                                  * the difference is the cumulative
 401                                  * total over the window size
 402                                  */
 403                                 temp_err = F_VAL(series, F_FIRST(series)) -
 404                                        F_VAL(series,windex);
 405                                 sq_err = temp_err / win;
 406                         }
 407                 }
 408                 else
 409                 {
 410                         sq_err = DBIG_VAL;
 411                 }
 412
 413                 if(sq_err < min_se)
 414                 {
 415                         min_i = i;
 416                         min_se = sq_err;
 417                 }
 418         }
 419
 420         /*
 421          * report the best
 422          */
 423         ferr =
 424 (l_forcl->forcs[min_i])->forecast((l_forcl->forcs[min_i])->state,&forecast);
 425
 426         /*
 427          * if an error has occurred, print and error and report
 428          * FORE_ERROR_VAL
 429          */
 430         if(ferr == 0)
 431         {
 432                 forecast = F_VAL(l_forcl->forcs[0]->series,
 433                                 F_FIRST(l_forcl->forcs[0]->series));
 434
 435                 return(0);
 436         }
 437
 438         /*
 439          * update the method that was used
 440          */
 441         s->method = min_i;
 442         *out_f = forecast;
 443
 444         return(1);
 445 }
 446
 447 int
 448 LocalWinMAEForecast(char *state, double *out_f)
 449 {
 450         struct mse_state *s;
 451         int i;
 452         int win;
 453         int min_i;
 454         double min_ae;
 455         double err;
 456         double forecast;
 457         forcl l_forcl;
 458         int ferr;
 459         double temp_err;
 460         int windex;
 461         fbuff series;
 462
 463         s = (struct mse_state *)state;
 464
 465         l_forcl = s->forc_list;
 466         win = s->win;
 467
 468
 469         min_ae = DBIG_VAL;
 470         min_i = 0;
 471
 472         /*
 473          * find the forecaster that has the lowest current
 474          * mean square error --
 475          *
 476          * if the win size > 0, look only over the window
 477          * using the cumulative series.  Otherwise, use the total
 478          * cumulative value
 479          */
 480         for(i=0; i < l_forcl->count; i++)
 481         {
 482                 if((l_forcl->forcs[i])->count > 0.0)
 483                 {
 484                         series =
 485                         (l_forcl->forcs[i])->ae_series;
 486
 487                         /*
 488                          * make sure we only window as much
 489                          * as there is
 490                          */
 491                         if(win > (F_COUNT(series) - 1))
 492                                 win = F_COUNT(series) - 1;
 493
 494                         if(win == 0)
 495                         {
 496                                 /*
 497                                  * get the first value which is
 498                                  * the current cumulative total
 499                                  */
 500                                 temp_err =
 501                                 F_VAL(series, F_FIRST(series));
 502                                 err = temp_err /
 503                                         (l_forcl->forcs[i])->count;
 504                         }
 505                         else
 506                         {
 507                                 /*
 508                                  * get the index of the end of the
 509                                  * window
 510                                  */
 511                                 windex = MODPLUS(F_FIRST(series),
 512                                                  win,
 513                                                  F_SIZE(series));
 514                                 /*
 515                                  * the difference is the cumulative
 516                                  * total over the window size
 517                                  */
 518                                 temp_err =
 519                                         F_VAL(series, F_FIRST(series)) -
 520                                         F_VAL(series,windex);
 521                                 err = temp_err / win;
 522                         }
 523                 }
 524                 else
 525                 {
 526                         err = DBIG_VAL;
 527                 }
 528
 529                 if(err < min_ae)
 530                 {
 531                         min_i = i;
 532                         min_ae = err;
 533                 }
 534         }
 535
 536         /*
 537          * report the best
 538          */
 539         ferr =
 540 (l_forcl->forcs[min_i])->forecast((l_forcl->forcs[min_i])->state,&forecast);
 541
 542         /*
 543          * if an error has occurred, print and error and report
 544          * FORE_ERROR_VAL
 545          */
 546         if(ferr == 0)
 547         {
 548                 forecast = F_VAL(l_forcl->forcs[0]->series,
 549                                 F_FIRST(l_forcl->forcs[0]->series));
 550
 551                 *out_f = forecast;
 552                 return(0);
 553         }
 554
 555         /*
 556          * update the method that was used
 557          */
 558         s->method = min_i;
 559         *out_f = forecast;
 560
 561         return(1);
 562 }
 563
 564 int
 565 DerivedMethod(char *state)
 566 {
 567         struct mse_state *s;
 568
 569         s = (struct mse_state *)state;
 570
 571         return(s->method);
 572 }
 573
 574 /*
 575  * lifetime versions
 576  */
 577 double
 578 MSELifetimeForecast(char *state)
 579 {
 580         forclife flife;
 581         forcl l_forcl;
 582
 583         flife = (forclife)state;
 584         l_forcl = (forcl)flife->forc_list;
 585         if(F_COUNT(l_forcl->forcs[0]->series) == 0)
 586         {
 587                 if(flife->count == 0.0)
 588                         return(0.0);
 589                 return(flife->total / flife->count);
 590         }
 591         return(MSEForecast((char *)flife->forc_list));
 592 }
 593
 594 double
 595 MSELifetimeError(char *state)
 596 {
 597         forclife flife;
 598         flife = (forclife)state;
 599         return(MSEError((char *)flife->forc_list));
 600 }
 601
 602 int
 603 MSELifetimeMethod(char *state)
 604 {
 605         forclife flife;
 606         flife = (forclife)state;
 607         return(MSEMethod((char *)flife->forc_list));
 608 }
 609
 610 double
 611 MSELifetimeEpoch(char *state)
 612 {
 613         forclife flife;
 614         flife = (forclife)state;
 615         return(flife->epoch_end);
 616 }
 617
 618 double
 619 MAELifetimeForecast(char *state)
 620 {
 621         forclife flife;
 622         forcl l_forcl;
 623
 624         flife = (forclife)state;
 625         l_forcl = (forcl)flife->forc_list;
 626
 627         if(F_COUNT(l_forcl->forcs[0]->series) == 0)
 628         {
 629                 if(flife->count == 0.0)
 630                         return(0.0);
 631                 return(flife->total / flife->count);
 632         }
 633         return(MAEForecast((char *)flife->forc_list));
 634 }
 635
 636 double
 637 MAELifetimeError(char *state)
 638 {
 639         forclife flife;
 640         flife = (forclife)state;
 641         return(MAEError((char *)flife->forc_list));
 642 }
 643
 644 int
 645 MAELifetimeMethod(char *state)
 646 {
 647         forclife flife;
 648         flife = (forclife)state;
 649         return(MAEMethod((char *)flife->forc_list));
 650 }
 651
 652 double
 653 MAELifetimeEpoch(char *state)
 654 {
 655         forclife flife;
 656         flife = (forclife)state;
 657         return(flife->epoch_end);
 658 }