src/nws_portability/Forecast/forc.c

   1 #include <unistd.h>
   2 #include <stdlib.h>
   3 #include <stdio.h>
   4 #include <math.h>
   5 #include <string.h>
   6 #include <strings.h>
   7
   8
   9
  10 /*
  11  * individual forecaster packages
  12  */
  13 #include "forecasters.h"
  14
  15 #include "mse_forc.h"
  16 #include "forc.h"
  17
  18 #define FORECASTER(name,init,params,freer,update,forecaster)    \
  19 {                                                               \
  20         frb = InitForcB(series,                                 \
  21                         time_stamps,                            \
  22                         (name),                                 \
  23                         (init),                                 \
  24                         (params),                               \
  25                         (freer),                                \
  26                         (update),                               \
  27                         (forecaster));                          \
  28                                                                 \
  29         if(frb == NULL)                                         \
  30         {                                                       \
  31                 free(l_forcl);                                  \
  32                 free(forcs);                                    \
  33                 free(derived_forcs);                            \
  34                 FreeFBuff(series);                              \
  35                 FreeFBuff(time_stamps);                         \
  36                 return(NULL);                                   \
  37         }                                                       \
  38                                                                 \
  39         frb->series = series;                                   \
  40         frb->time_stamps = time_stamps;                         \
  41         forcs[l_forcl->count] = frb;                            \
  42         l_forcl->count = l_forcl->count + 1;                    \
  43 }
  44
  45
  46 #define DERIVED_FORC(name,init,params,freer,update,forecaster)  \
  47 {                                                               \
  48         frb = InitForcB(series,                                 \
  49                   time_stamps,                                  \
  50                   (name),                                       \
  51                   (init),                                       \
  52                   (params),                                     \
  53                   (freer),                                      \
  54                   (update),                                     \
  55                   (forecaster));                                \
  56                                                                 \
  57         if(frb == NULL)                                         \
  58         {                                                       \
  59                 free((char *)forcs);                                    \
  60                 free((char *)derived_forcs);                            \
  61                 FreeFBuff(series);                              \
  62                 FreeFBuff(time_stamps);                         \
  63                 for(i=0; i < l_forcl->count; i++)               \
  64                 {                                               \
  65                         FreeForcB(l_forcl->forcs[i]);           \
  66                 }                                               \
  67                 free((char *)l_forcl);                                  \
  68                 return(NULL);                                   \
  69         }                                                       \
  70         frb->series = series;                                   \
  71         frb->time_stamps = time_stamps;                         \
  72         derived_forcs[l_forcl->derived_count] = frb;            \
  73         l_forcl->derived_count = l_forcl->derived_count + 1;    \
  74 }
  75
  76 forcb
  77 InitForcB(fbuff series,
  78           fbuff time_stamps,
  79           const char *name,
  80           char *(*init)(),
  81           const char *params,
  82           void (*freer)(),
  83           void (*update)(),
  84           int (*forecast)())
  85 {
  86         forcb frb;
  87
  88         frb = (forcb)malloc(FORCB_SIZE);
  89
  90
  91         if(frb == NULL)
  92         {
  93                 return(NULL);
  94         }
  95
  96         /*
  97          * get fbuff space for the cumulative error series -- needed
  98          * for windowed MSE and MAE predictions
  99          */
 100         frb->se_series = InitFBuff(MAX_MSE_WIN);
 101
 102         if(frb->se_series == NULL)
 103         {
 104                 free(frb);
 105                 return(NULL);
 106         }
 107
 108         frb->ae_series = InitFBuff(MAX_MSE_WIN);
 109
 110         if(frb->se_series == NULL)
 111         {
 112                 FreeFBuff(frb->se_series);
 113                 free(frb);
 114                 return(NULL);
 115         }
 116         /*
 117          * could be NULL because of derived forecasters
 118          */
 119         if(init != NULL)
 120         {
 121                 frb->state = (*init)(series,time_stamps,params);
 122                 if(frb == NULL)
 123                 {
 124                         FreeFBuff(frb->se_series);
 125                         FreeFBuff(frb->ae_series);
 126                         free(frb);
 127                         fprintf(stderr,"InitForcB: init failed for %s\n",
 128                                         name);
 129                         fflush(stderr);
 130                         return(NULL);
 131                 }
 132         }
 133         else
 134         {
 135                 frb->state = NULL;
 136         }
 137
 138         frb->update = update;
 139         frb->forecast = forecast;
 140         frb->free = freer;
 141
 142
 143         frb->se = 0.0;
 144         frb->ae = 0.0;
 145         frb->count = 0.0;
 146         strncpy(frb->name,name,sizeof(frb->name));
 147
 148         /*
 149          * point forcb at the series and time stamp buffers
 150          */
 151         frb->series = series;
 152         frb->time_stamps = time_stamps;
 153
 154         frb->best_f = 0.0;
 155         frb->best_i = 0;
 156         frb->best_err = 0.0;
 157
 158
 159         return(frb);
 160 }
 161
 162 void
 163 FreeForcB(forcb frb)
 164 {
 165         if(frb->state != NULL)
 166         {
 167                 frb->free(frb->state);
 168         }
 169
 170         /*
 171          * defensive programming
 172          */
 173         if(frb->se_series != NULL)
 174         {
 175                 FreeFBuff(frb->se_series);
 176         }
 177         if(frb->ae_series != NULL)
 178         {
 179                 FreeFBuff(frb->ae_series);
 180         }
 181         free(frb);
 182
 183         return;
 184 }
 185
 186 void
 187 FreeForcl(char *i_forcl)
 188 {
 189         int i;
 190         forcl l_forcl;
 191
 192         l_forcl = (forcl)i_forcl;
 193
 194         /*
 195          * shared by all forcbs
 196          */
 197         FreeFBuff(l_forcl->forcs[0]->series);
 198         FreeFBuff(l_forcl->forcs[0]->time_stamps);
 199
 200         for(i=0; i < l_forcl->count; i++)
 201         {
 202                 FreeForcB(l_forcl->forcs[i]);
 203         }
 204         for(i=0; i < l_forcl->derived_count; i++)
 205         {
 206                 FreeForcB(l_forcl->derived_forcs[i]);
 207         }
 208
 209         /*
 210          * the state is the l_forcl itself -- clear it to free works
 211          * properly
 212          */
 213         l_forcl->total_mse->state = NULL;
 214         l_forcl->total_mae->state = NULL;
 215         FreeForcB(l_forcl->total_mse);
 216         FreeForcB(l_forcl->total_mae);
 217
 218         free(l_forcl->derived_forcs);
 219         free(l_forcl->forcs);
 220         free(l_forcl);
 221
 222         return;
 223 }
 224
 225 /*
 226  * initializes a single fbuff for the data series and another for the
 227  * time stamps.
 228  *
 229  * initializes an array of forcb structs, one per forecaster
 230  *
 231  * points each forcb to the series and time stamps
 232  */
 233 char *
 234 InitForcl(int max_forc_count, int buff_size)
 235 {
 236         forcb frb;
 237         forcb *forcs;
 238         forcb *derived_forcs;
 239         forcl l_forcl;
 240         fbuff series;
 241         fbuff time_stamps;
 242         int i;
 243         char derived_params[255];
 244
 245         /*
 246          * first, get space for series and time stamps
 247          */
 248         series = InitFBuff(buff_size);
 249         if(series == NULL)
 250         {
 251                 fprintf(stderr,"InitForcs: couldn't get fbuff for series\n");
 252                 fflush(stderr);
 253                 return(NULL);
 254         }
 255
 256         time_stamps = InitFBuff(buff_size);
 257         if(time_stamps == NULL)
 258         {
 259                 fprintf(stderr,
 260                         "InitForcs: couldn't get fbuff for time stamps\n");
 261                 fflush(stderr);
 262                 FreeFBuff(series);
 263                 return(NULL);
 264         }
 265
 266         /*
 267          * now, get forcl list space
 268          */
 269         l_forcl = (forcl)(malloc(FORCL_SIZE));
 270         if(l_forcl == NULL)
 271         {
 272                 fprintf(stderr,"InitForcs: couldn't malloc %d forcl\n",
 273                                 max_forc_count);
 274                 fflush(stderr);
 275                 FreeFBuff(series);
 276                 FreeFBuff(time_stamps);
 277                 return(NULL);
 278         }
 279
 280         memset((char *)l_forcl,0,FORCL_SIZE);
 281         l_forcl->derived_count = 0;
 282
 283         /*
 284          * get space for individual forecasters
 285          */
 286         forcs = (forcb *)(malloc(max_forc_count*sizeof(forcb)));
 287
 288         if(forcs == NULL)
 289         {
 290                 fprintf(stderr,"InitForcl: couldn't malloc %d forcb pointers\n",
 291                                 max_forc_count);
 292                 fflush(stderr);
 293                 free(l_forcl);
 294                 FreeFBuff(series);
 295                 FreeFBuff(time_stamps);
 296                 return(NULL);
 297         }
 298
 299         /*
 300          * get space for derived forecasters
 301          */
 302         derived_forcs = (forcb *)(malloc(max_forc_count*sizeof(forcb)));
 303
 304         if(derived_forcs == NULL)
 305         {
 306                 fprintf(stderr,"InitForcl: couldn't malloc %d derived forcb pointers\n",
 307                                 max_forc_count);
 308                 fflush(stderr);
 309                 free(l_forcl);
 310                 FreeFBuff(series);
 311                 FreeFBuff(time_stamps);
 312                 free(forcs);
 313                 return(NULL);
 314         }
 315
 316         l_forcl->forcs = forcs;
 317         l_forcl->derived_forcs = derived_forcs;
 318
 319         /*
 320          * include the primary forecasters
 321          */
 322 #include "forecasters.c"
 323
 324
 325         /*
 326          * now do derived forecaster types
 327          */
 328
 329         /*
 330          * for min MSE, we set the window size to 0.  The first
 331          * argument is the address of the forc list structure we are building
 332          *
 333          * thread safe?  I think not.
 334          */
 335         sprintf(derived_params,
 336                 "%p %d",
 337                 (void *)l_forcl,
 338                 0);
 339         DERIVED_FORC("Minimum MSE",
 340                      InitWinMSE,
 341                      derived_params,
 342                      FreeWinMSE,
 343                      NULL,
 344                      LocalWinMSEForecast);
 345
 346         /*
 347          * for min MAE, we set the window size to 0.  The first
 348          * argument is the address of the forc list structure we are building
 349          *
 350          * thread safe?  I think not.
 351          */
 352         sprintf(derived_params,
 353                 "%p %d",
 354                 (void *)l_forcl,
 355                 0);
 356
 357         DERIVED_FORC("Minimum MAE",
 358                      InitWinMAE,
 359                      derived_params,
 360                      FreeWinMAE,
 361                      NULL,
 362                      LocalWinMAEForecast);
 363
 364         /*
 365          * for min windowed MSE, we set the window size to 1.  The first
 366          * argument is the address of the forc list structure we are building
 367          *
 368          * thread safe?  I think not.
 369          */
 370         sprintf(derived_params,
 371                 "%p %d",
 372                 (void *)l_forcl,
 373                 1);
 374         DERIVED_FORC("Minimum Window 1 MSE",
 375                      InitWinMSE,
 376                      derived_params,
 377                      FreeWinMSE,
 378                      NULL,
 379                      LocalWinMSEForecast);
 380         /*
 381          * for min windowed MSE, we set the window size to 5.  The first
 382          * argument is the address of the forc list structure we are building
 383          *
 384          * thread safe?  I think not.
 385          */
 386         sprintf(derived_params,
 387                 "%p %d",
 388                 (void *)l_forcl,
 389                 5);
 390         DERIVED_FORC("Minimum Window 5 MSE",
 391                      InitWinMSE,
 392                      derived_params,
 393                      FreeWinMSE,
 394                      NULL,
 395                      LocalWinMSEForecast);
 396
 397         /*
 398          * for min windowed MSE, we set the window size to 10.  The first
 399          * argument is the address of the forc list structure we are building
 400          *
 401          * thread safe?  I think not.
 402          */
 403         sprintf(derived_params,
 404                 "%p %d",
 405                 (void *)l_forcl,
 406                 10);
 407         DERIVED_FORC("Minimum Window 10 MSE",
 408                      InitWinMSE,
 409                      derived_params,
 410                      FreeWinMSE,
 411                      NULL,
 412                      LocalWinMSEForecast);
 413
 414         /*
 415          * for min windowed MSE, we set the window size to 20.  The first
 416          * argument is the address of the forc list structure we are building
 417          *
 418          * thread safe?  I think not.
 419          */
 420         sprintf(derived_params,
 421                 "%p %d",
 422                 (void *)l_forcl,
 423                 20);
 424         DERIVED_FORC("Minimum Window 20 MSE",
 425                      InitWinMSE,
 426                      derived_params,
 427                      FreeWinMSE,
 428                      NULL,
 429                      LocalWinMSEForecast);
 430         /*
 431          * for min windowed MSE, we set the window size to 30.  The first
 432          * argument is the address of the forc list structure we are building
 433          *
 434          * thread safe?  I think not.
 435          */
 436         sprintf(derived_params,
 437                 "%p %d",
 438                 (void *)l_forcl,
 439                 30);
 440         DERIVED_FORC("Minimum Window 30 MSE",
 441                      InitWinMSE,
 442                      derived_params,
 443                      FreeWinMSE,
 444                      NULL,
 445                      LocalWinMSEForecast);
 446         /*
 447          * for min windowed MSE, we set the window size to 50.  The first
 448          * argument is the address of the forc list structure we are building
 449          *
 450          * thread safe?  I think not.
 451          */
 452         sprintf(derived_params,
 453                 "%p %d",
 454                 (void *)l_forcl,
 455                 50);
 456         DERIVED_FORC("Minimum Window 50 MSE",
 457                      InitWinMSE,
 458                      derived_params,
 459                      FreeWinMSE,
 460                      NULL,
 461                      LocalWinMSEForecast);
 462         /*
 463          * for min windowed MSE, we set the window size to 100.  The first
 464          * argument is the address of the forc list structure we are building
 465          *
 466          * thread safe?  I think not.
 467          */
 468         sprintf(derived_params,
 469                 "%p %d",
 470                 (void *)l_forcl,
 471                 100);
 472         DERIVED_FORC("Minimum Window 100 MSE",
 473                      InitWinMSE,
 474                      derived_params,
 475                      FreeWinMSE,
 476                      NULL,
 477                      LocalWinMSEForecast);
 478
 479         /*
 480          * for win min MAE, with window size to 1.  The first
 481          * argument is the address of the forc list structure we are building
 482          *
 483          * thread safe?  I think not.
 484          */
 485         sprintf(derived_params,
 486                 "%p %d",
 487                 (void *)l_forcl,
 488                 1);
 489
 490         DERIVED_FORC("Minimum Window 1 MAE",
 491                      InitWinMAE,
 492                      derived_params,
 493                      FreeWinMAE,
 494                      NULL,
 495                      LocalWinMAEForecast);
 496         /*
 497          * for win min MAE, with window size to 5.  The first
 498          * argument is the address of the forc list structure we are building
 499          *
 500          * thread safe?  I think not.
 501          */
 502         sprintf(derived_params,
 503                 "%p %d",
 504                 (void *)l_forcl,
 505                 5);
 506
 507         DERIVED_FORC("Minimum Window 5 MAE",
 508                      InitWinMAE,
 509                      derived_params,
 510                      FreeWinMAE,
 511                      NULL,
 512                      LocalWinMAEForecast);
 513         /*
 514          * for win min MAE, with window size to 10.  The first
 515          * argument is the address of the forc list structure we are building
 516          *
 517          * thread safe?  I think not.
 518          */
 519         sprintf(derived_params,
 520                 "%p %d",
 521                 (void *)l_forcl,
 522                 10);
 523
 524         DERIVED_FORC("Minimum Window 10 MAE",
 525                      InitWinMAE,
 526                      derived_params,
 527                      FreeWinMAE,
 528                      NULL,
 529                      LocalWinMAEForecast);
 530         /*
 531          * for win min MAE, with window size to 20.  The first
 532          * argument is the address of the forc list structure we are building
 533          *
 534          * thread safe?  I think not.
 535          */
 536         sprintf(derived_params,
 537                 "%p %d",
 538                 (void *)l_forcl,
 539                 20);
 540
 541         DERIVED_FORC("Minimum Window 20 MAE",
 542                      InitWinMAE,
 543                      derived_params,
 544                      FreeWinMAE,
 545                      NULL,
 546                      LocalWinMAEForecast);
 547         /*
 548          * for win min MAE, with window size to 30.  The first
 549          * argument is the address of the forc list structure we are building
 550          *
 551          * thread safe?  I think not.
 552          */
 553         sprintf(derived_params,
 554                 "%p %d",
 555                 (void *)l_forcl,
 556                 30);
 557
 558         DERIVED_FORC("Minimum Window 30 MAE",
 559                      InitWinMAE,
 560                      derived_params,
 561                      FreeWinMAE,
 562                      NULL,
 563                      LocalWinMAEForecast);
 564         /*
 565          * for win min MAE, with window size to 50.  The first
 566          * argument is the address of the forc list structure we are building
 567          *
 568          * thread safe?  I think not.
 569          */
 570         sprintf(derived_params,
 571                 "%p %d",
 572                 (void *)l_forcl,
 573                 50);
 574
 575         DERIVED_FORC("Minimum Window 50 MAE",
 576                      InitWinMAE,
 577                      derived_params,
 578                      FreeWinMAE,
 579                      NULL,
 580                      LocalWinMAEForecast);
 581         /*
 582          * for win min MAE, with window size to 100.  The first
 583          * argument is the address of the forc list structure we are building
 584          *
 585          * thread safe?  I think not.
 586          */
 587         sprintf(derived_params,
 588                 "%p %d",
 589                 (void *)l_forcl,
 590                 100);
 591
 592         DERIVED_FORC("Minimum Window 100 MAE",
 593                      InitWinMAE,
 594                      derived_params,
 595                      FreeWinMAE,
 596                      NULL,
 597                      LocalWinMAEForecast);
 598         /*
 599          * end of derived forecaster init
 600          */
 601
 602         /*
 603          * now init a total MSE forcb for the
 604          * list
 605          */
 606
 607         frb = InitForcB(series,
 608                         time_stamps,
 609                         "Total MSE",
 610                         NULL,           /* no init */
 611                         NULL,           /* no params */
 612                         NULL,           /* no free */
 613                         NULL,           /* no update */
 614                         TotalMSEForecast);      /* MSE forecast */
 615         if(frb == NULL)
 616         {
 617                 free(forcs);
 618                 free(derived_forcs);
 619                 FreeFBuff(series);
 620                 FreeFBuff(time_stamps);
 621                 for(i=0; i < l_forcl->count; i++)
 622                 {
 623                         FreeForcB(l_forcl->forcs[i]);
 624                 }
 625                 for(i=0; i < l_forcl->derived_count; i++)
 626                 {
 627                         FreeForcB(l_forcl->derived_forcs[i]);
 628                 }
 629                 free(l_forcl->derived_forcs);
 630                 free(l_forcl);
 631                 return(NULL);
 632         }
 633         frb->series = series;
 634         frb->time_stamps = time_stamps;
 635         frb->ae = 0.0;
 636         frb->se = 0.0;
 637         frb->count = 0.0;
 638         /*
 639          * total forecasters take forcl as state
 640          */
 641         frb->state = (char *)l_forcl;
 642         l_forcl->total_mse = frb;
 643
 644         /*
 645          * now init a total MAE forcb for the
 646          * list
 647          */
 648
 649         frb = InitForcB(series,
 650                         time_stamps,
 651                         "Total MAE",
 652                         NULL,           /* no init */
 653                         NULL,           /* no params */
 654                         NULL,           /* no free */
 655                         NULL,           /* no update */
 656                         TotalMAEForecast);      /* MSE forecast */
 657         if(frb == NULL)
 658         {
 659                 free(forcs);
 660                 free(derived_forcs);
 661                 FreeFBuff(series);
 662                 FreeFBuff(time_stamps);
 663                 for(i=0; i < l_forcl->count; i++)
 664                 {
 665                         FreeForcB(l_forcl->forcs[i]);
 666                 }
 667                 for(i=0; i < l_forcl->derived_count; i++)
 668                 {
 669                         FreeForcB(l_forcl->derived_forcs[i]);
 670                 }
 671                 free(l_forcl->derived_forcs);
 672                 free(l_forcl);
 673                 return(NULL);
 674         }
 675         frb->series = series;
 676         frb->time_stamps = time_stamps;
 677         frb->ae = 0.0;
 678         frb->se = 0.0;
 679         frb->count = 0.0;
 680         /*
 681          * total forecasters take forcl as state
 682          */
 683         frb->state = (char *)l_forcl;
 684         l_forcl->total_mae = frb;
 685
 686
 687
 688         return((char *)l_forcl);
 689
 690 }
 691
 692 void
 693 UpdateForecasts(char *i_forcl, double ts, double value)
 694 {
 695         int i;
 696         double err;
 697         double forecast;
 698         int ferr;
 699         forcl l_forcl;
 700         int forc_okay;
 701         double best_sq_f;
 702         double best_f;
 703         int best_i;
 704         int best_sq_i;
 705         double best_sq_err;
 706         double best_err;
 707
 708         l_forcl = (forcl)i_forcl;
 709
 710         forc_okay = 0;
 711
 712         /*
 713          * get mse forecast for this value and update error
 714          */
 715         ferr = (l_forcl->total_mse)->forecast((l_forcl->total_mse)->state,
 716                               &forecast);
 717         if(ferr == 1)
 718         {
 719                 forc_okay = 1;
 720                 err = forecast - value;
 721                 if(err < 0.0)
 722                 {
 723                         err = err * -1.0;
 724                 }
 725                 l_forcl->total_mse->ae += err;
 726                 l_forcl->total_mse->se += err*err;
 727                 l_forcl->total_mse->count += 1.0;
 728
 729                 /*
 730                  * write out the cumulative total in the
 731                  * series buffer so that we may do windowed
 732                  * mse and mae things
 733                  */
 734                 UpdateFBuff(l_forcl->total_mse->se_series,
 735                            l_forcl->total_mse->se);
 736                 UpdateFBuff(l_forcl->total_mse->ae_series,
 737                            l_forcl->total_mse->ae);
 738         }
 739         /*
 740          * get mae forecast for this value and update error
 741          */
 742         ferr = (l_forcl->total_mae)->forecast((l_forcl->total_mae)->state,
 743                               &forecast);
 744         if(ferr == 1)
 745         {
 746                 err = forecast - value;
 747                 if(err < 0.0)
 748                 {
 749                         err = err * -1.0;
 750                 }
 751                 l_forcl->total_mae->ae += err;
 752                 l_forcl->total_mae->se += err*err;
 753                 l_forcl->total_mae->count += 1.0;
 754
 755                 /*
 756                  * write out the cumulative total in the
 757                  * series buffer so that we may do windowed
 758                  * mse and mae things
 759                  */
 760                 UpdateFBuff(l_forcl->total_mae->se_series,
 761                            l_forcl->total_mae->se);
 762                 UpdateFBuff(l_forcl->total_mae->ae_series,
 763                            l_forcl->total_mae->ae);
 764         }
 765
 766         if(forc_okay == 1)
 767         {
 768                 /*
 769                  * update derived forecast error values
 770                  */
 771                 for(i=0; i < l_forcl->derived_count; i++)
 772                 {
 773                         /*
 774                          * get forecast for this value and update error
 775                          */
 776                         ferr =
 777         (l_forcl->derived_forcs[i])->forecast((l_forcl->derived_forcs[i])->state,
 778                                               &forecast);
 779                         if(ferr == 0)
 780                         {
 781                                 continue;
 782                         }
 783                         forc_okay = 1;
 784                         err = forecast - value;
 785                         if(err < 0.0)
 786                         {
 787                                 err = err * -1.0;
 788                         }
 789                         l_forcl->derived_forcs[i]->ae += err;
 790                         l_forcl->derived_forcs[i]->se += err*err;
 791                         l_forcl->derived_forcs[i]->count += 1.0;
 792
 793                         /*
 794                          * write out the cumulative total in the
 795                          * series buffer so that we may do windowed
 796                          * mse and mae things
 797                          */
 798                         UpdateFBuff(l_forcl->derived_forcs[i]->se_series,
 799                                    l_forcl->derived_forcs[i]->se);
 800                         UpdateFBuff(l_forcl->derived_forcs[i]->ae_series,
 801                                    l_forcl->derived_forcs[i]->ae);
 802                 }
 803         }
 804         /*
 805          * first, update the error values
 806          */
 807         best_sq_err = DBIG_VAL;
 808         best_err = DBIG_VAL;
 809         best_sq_f = 0.0;
 810         best_f = 0.0;
 811         best_sq_i = 0;
 812         best_i = 0;
 813         for(i=0; i < l_forcl->count; i++)
 814         {
 815                 /*
 816                  * get forecast for this value and update error
 817                  */
 818                 ferr =
 819         (l_forcl->forcs[i])->forecast((l_forcl->forcs[i])->state,&forecast);
 820                 if(ferr == 0)
 821                 {
 822                         continue;
 823                 }
 824
 825                 err = forecast - value;
 826                 if(err < 0.0)
 827                 {
 828                         err = err * -1.0;
 829                 }
 830                 l_forcl->forcs[i]->ae += err;
 831                 l_forcl->forcs[i]->se += err*err;
 832                 l_forcl->forcs[i]->count += 1.0;
 833
 834                 /*
 835                  * find the forecast that actually has the min error
 836                  */
 837                 if((err*err) < best_sq_err)
 838                 {
 839                         best_sq_err = err*err;
 840                         best_sq_i = i;
 841                         best_sq_f = forecast;
 842                 }
 843                 if(err < best_err)
 844                 {
 845                         best_err = err;
 846                         best_i = i;
 847                         best_f = forecast;
 848                 }
 849                 /*
 850                  * write out the cumulative total in the
 851                  * series buffer so that we may do windowed
 852                  * mse and mae things
 853                  */
 854                 UpdateFBuff(l_forcl->forcs[i]->se_series,
 855                            l_forcl->forcs[i]->se);
 856                 UpdateFBuff(l_forcl->forcs[i]->ae_series,
 857                            l_forcl->forcs[i]->ae);
 858         }
 859
 860         /*
 861          * record the forecast that actually won
 862          */
 863         if(best_sq_err != DBIG_VAL)
 864         {
 865                 l_forcl->total_mse->best_f = best_sq_f;
 866                 l_forcl->total_mse->best_i = best_sq_i;
 867                 l_forcl->total_mse->best_err += best_sq_err;
 868
 869                 l_forcl->total_mae->best_f = best_f;
 870                 l_forcl->total_mae->best_i = best_i;
 871                 l_forcl->total_mae->best_err += best_err;
 872         }
 873
 874         /*
 875          * now update the series and time_stamps fbuffs.  all forcs point
 876          * the same series and time_stamp fbuffs
 877          */
 878
 879         if(l_forcl->count >= 1)
 880         {
 881                 UpdateFBuff(l_forcl->forcs[0]->series, value);
 882                 UpdateFBuff(l_forcl->forcs[0]->time_stamps, ts);
 883         }
 884
 885         for(i=0; i < l_forcl->count; i++)
 886         {
 887                 /*
 888                  * now update forecast state for the individual forecasters
 889                  *
 890                  * note that the forecasters need not record the series and
 891                  * time_stamp points, hence the ts and value parameters
 892                  */
 893                 (l_forcl->forcs[i])->update((l_forcl->forcs[i])->state,
 894                                    ts,
 895                                    value);
 896         }
 897
 898         return;
 899 }
 900
 901 int
 902 ForcRange(char *i_forcl, double *low, double *high, int *low_i, int *high_i)
 903 {
 904         forcl l_forcl;
 905         int i;
 906         double min = DBIG_VAL;
 907         double max = -1.0 * DBIG_VAL;
 908         double f;
 909         int err;
 910
 911         l_forcl = (forcl)i_forcl;
 912
 913         /*
 914          * simple forecasters
 915          */
 916         for(i=0; i < l_forcl->count; i++)
 917         {
 918                 err = l_forcl->forcs[i]->forecast(l_forcl->forcs[i]->state,
 919                                                   &f);
 920                 if(err == 0)
 921                         continue;
 922
 923                 if(f < min)
 924                 {
 925                         min = f;
 926                         *low_i = i;
 927                 }
 928
 929                 if(f > max)
 930                 {
 931                         max = f;
 932                         *high_i = i;
 933                 }
 934         }
 935
 936         if((min == DBIG_VAL) || (max == (-1.0*DBIG_VAL)))
 937         {
 938                 return(0);
 939         }
 940
 941         *low = min;
 942         *high = max;
 943
 944         return(1);
 945
 946 }
 947
 948 void
 949 PrintForecastSummary(char *i_forcl)
 950 {
 951         forcl l_forcl;
 952         int i;
 953
 954         l_forcl = (forcl)i_forcl;
 955
 956         /*
 957          * simple forecasters
 958          */
 959         for(i=0; i < l_forcl->count; i++)
 960         {
 961                 fprintf(stdout,"MSE: %3.4f\tMAE: %3.4f\t%s\n",
 962                         l_forcl->forcs[i]->se/l_forcl->forcs[i]->count,
 963                         l_forcl->forcs[i]->ae/l_forcl->forcs[i]->count,
 964                         l_forcl->forcs[i]->name);
 965         }
 966
 967         /*
 968          * secondary forecasters
 969          */
 970         for(i=0; i < l_forcl->derived_count; i++)
 971         {
 972                 fprintf(stdout,"MSE: %3.4f\tMAE: %3.4f\t%s\n",
 973                 l_forcl->derived_forcs[i]->se/l_forcl->derived_forcs[i]->count,
 974                 l_forcl->derived_forcs[i]->ae/l_forcl->derived_forcs[i]->count,
 975                 l_forcl->derived_forcs[i]->name);
 976         }
 977
 978         /*
 979          * totals
 980          */
 981         fprintf(stdout,"MSE: %3.4f\tMAE: %3.4f\tTotal MSE\n",
 982                         l_forcl->total_mse->se / l_forcl->total_mse->count,
 983                         l_forcl->total_mse->ae / l_forcl->total_mse->count);
 984         fprintf(stdout,"MSE: %3.4f\tMAE: %3.4f\tTotal MAE\n",
 985                         l_forcl->total_mae->se / l_forcl->total_mae->count,
 986                         l_forcl->total_mae->ae / l_forcl->total_mae->count);
 987         fprintf(stdout,"MSE: %3.4f\tMAE: %3.4f\tOptimum\n",
 988                         l_forcl->total_mse->best_err/l_forcl->total_mse->count,
 989                         l_forcl->total_mae->best_err/l_forcl->total_mse->count);
 990
 991         return;
 992 }
 993
 994 void
 995 PrintLifetimeForecastSummary(char *state)
 996 {
 997         forclife flife;
 998
 999         flife = (forclife)state;
1000
1001         PrintForecastSummary(flife->forc_list);
1002
1003         return;
1004 }
1005
1006 /*
1007  * routine put in to support FORECASTAPI_MethodName which doesn't
1008  * take a state record
1009  */
1010 void
1011 GetForcNames(char *state, char *methodNames[], int max_size, int *out_size)
1012 {
1013         forcl s;
1014         int i;
1015
1016         s = (forcl)state;
1017
1018         for(i=0; i < max_size; i++)
1019         {
1020                 if(i >= s->count)
1021                         break;
1022
1023                 methodNames[i] = s->forcs[i]->name;
1024         }
1025
1026         *out_size = i;
1027
1028         return;
1029 }
1030
1031 char *
1032 InitForcLife(int max_forc_count, int buff_size, double lifetime)
1033 {
1034         char *l_forcl;
1035         forclife flife;
1036
1037         /*
1038          * get a forc list for the series
1039          */
1040         l_forcl = InitForcl(max_forc_count, buff_size);
1041         if(l_forcl == NULL)
1042         {
1043                 fprintf(stderr,"InitForcLife: no space for forc list\n");
1044                 fflush(stderr);
1045                 return(NULL);
1046         }
1047
1048         /*
1049          * now, get space for the lifetime forecasting structure
1050          */
1051         flife = (forclife)malloc(FORCLIFE_SIZE);
1052         if(flife == NULL)
1053         {
1054                 FreeForcl(l_forcl);
1055                 fprintf(stderr,"InitForcLife: no space for forclife struct\n");
1056                 fflush(stderr);
1057                 return(NULL);
1058         }
1059
1060
1061         flife->lifetime = lifetime;
1062         flife->epoch_end = 0.0;
1063         flife->total = 0.0;
1064         flife->count = 0.0;
1065         flife->forc_list = l_forcl;
1066
1067         return((char *)flife);
1068 }
1069
1070 void
1071 FreeForcLife(char *state)
1072 {
1073         forclife flife;
1074
1075         flife = (forclife)state;
1076
1077         FreeForcl(flife->forc_list);
1078         free(flife);
1079
1080         return;
1081 }
1082
1083 void
1084 UpdateForcLife(char *state, double ts, double value)
1085 {
1086         forclife flife;
1087         double avg;
1088
1089         flife = (forclife)state;
1090
1091         /*
1092          * if we are initializing
1093          */
1094         if(flife->epoch_end == 0.0)
1095         {
1096                 flife->epoch_end = ts + flife->lifetime;
1097         }
1098
1099         if(ts < flife->epoch_end)
1100         {
1101                 flife->total += value;
1102                 flife->count += 1.0;
1103         }
1104         else
1105         {
1106                 if(flife->count == 0.0)
1107                         flife->count = 1.0;
1108
1109                 avg = flife->total / flife->count;
1110                 UpdateForecasts((char *)flife->forc_list,
1111                                   flife->epoch_end,
1112                                   avg);
1113                 flife->total = 0.0;
1114                 flife->count = 0.0;
1115                 flife->epoch_end = ts+flife->lifetime;
1116         }
1117
1118         return;
1119 }
1120
1121 double
1122 LifetimeValue(char *state)
1123 {
1124         forclife flife;
1125         forcl l_forcl;
1126         double value;
1127
1128         flife = (forclife)state;
1129         l_forcl = (forcl)flife->forc_list;
1130
1131         value =
1132            F_VAL(l_forcl->forcs[0]->series,F_FIRST(l_forcl->forcs[0]->series));
1133
1134         return(value);
1135 }
1136
1137 double
1138 LifetimeTimestamp(char *state)
1139 {
1140         forclife flife;
1141         forcl l_forcl;
1142         double ts;
1143
1144         flife = (forclife)state;
1145         l_forcl = (forcl)flife->forc_list;
1146
1147         ts =
1148            F_VAL(l_forcl->forcs[0]->time_stamps,
1149                            F_FIRST(l_forcl->forcs[0]->time_stamps));
1150
1151         return(ts);
1152 }