src/nws_portability/Forecast/median.c

   1 /* $Id$ */
   2
   3 #include <unistd.h>
   4 #include <stdlib.h>
   5 #include <sys/types.h>
   6 #include <sys/times.h>
   7 #include <math.h>
   8 #include <stdio.h>
   9 #include <string.h>
  10
  11
  12 #include "fbuff.h"
  13 #include "median.h"
  14
  15
  16 /*
  17  * to do --
  18  *
  19  * init M_ts to BIG_TS
  20  *
  21  * handle the case where data is filling the window
  22  *
  23  * find the biggest window
  24  */
  25
  26 #define SMALL_TS (-1)
  27 #define BIG_TS (99999999999999999.99)
  28
  29
  30
  31 /*
  32  * state record for median and trimmed median
  33  */
  34 struct median_state
  35 {
  36         fbuff series;
  37         fbuff time_stamps;
  38         double M_array[MAX_MED_SIZE];
  39         double M_ts[MAX_MED_SIZE];
  40         int artificial_time;
  41         int M_count;
  42         int M_size;
  43         double trim;
  44 };
  45
  46 /*
  47  * state record for adjusted median and adjusted mean
  48  */
  49 struct ad_median_state
  50 {
  51         fbuff series;
  52         fbuff time_stamps;
  53         double M_array[MAX_MED_SIZE];
  54         double M_ts[MAX_MED_SIZE];
  55         int M_count;
  56         int artificial_time;
  57         int last_win;
  58         int win;
  59         int min;
  60         int max;
  61         int offset;
  62 };
  63
  64 void
  65 MSort(double *M_value, double *M_ts, double value, double ts, int size)
  66 {
  67         int i;
  68         int prev;
  69         int curr;
  70         int next;
  71         int l_index = 0;
  72         double temp;
  73         double min_ts = BIG_TS;
  74
  75         /*
  76          * list is sorted from smallest to biggest -- find the
  77          * oldest time stamp
  78          */
  79         for(i=0; i < size; i++)
  80         {
  81                 if(M_ts[i] < min_ts)
  82                 {
  83                         l_index = i;
  84                         min_ts = M_ts[i];
  85                 }
  86         }
  87
  88         /*
  89          * overwrite the slot with the new value and
  90          * timestamp
  91          */
  92         M_value[l_index] = value;
  93         M_ts[l_index] = ts;
  94
  95         /*
  96          * bump the value into the right place in the list
  97          */
  98
  99         prev = l_index - 1;
 100         if(prev < 0)
 101                 prev = 0;
 102         curr = l_index;
 103         next = l_index + 1;
 104         if(next > (size - 1))
 105                 next = size - 1;
 106         /*
 107          * while the value is not in the right spot
 108          */
 109         while(!((M_value[curr] >= M_value[prev]) &&
 110                 (M_value[curr] <= M_value[next])))
 111         {
 112                 /*
 113                  * if it is too big, shift up
 114                  */
 115                 if(M_value[curr] < M_value[prev])
 116                 {
 117                         temp = M_value[prev];
 118                         M_value[prev] = M_value[curr];
 119                         M_value[curr] = temp;
 120
 121                         temp = M_ts[prev];
 122                         M_ts[prev] = M_ts[curr];
 123                         M_ts[curr] = temp;
 124
 125                         curr = curr - 1;
 126                 }
 127                 else /* shift down */
 128                 {
 129                         temp = M_value[next];
 130                         M_value[next] = M_value[curr];
 131                         M_value[curr] = temp;
 132
 133                         temp = M_ts[next];
 134                         M_ts[next] = M_ts[curr];
 135                         M_ts[curr] = temp;
 136
 137                         curr = curr + 1;
 138                 }
 139
 140                 prev = curr - 1;
 141                 if(prev < 0)
 142                         prev = 0;
 143
 144                 next = curr + 1;
 145                 if(next > (size - 1))
 146                         next = size - 1;
 147         }
 148
 149         /*
 150          * new list should be sorted
 151          */
 152
 153         return;
 154 }
 155
 156
 157 /*
 158  * these next two routines can do windowed medians for different
 159  * window sizes.  We put them in for the Adjusted versions but
 160  * all medians should use them
 161  */
 162 double
 163 FindMedian(double *sorted_list,
 164            double *sorted_ts,
 165            int list_size,
 166            double now,
 167            int win_size)
 168 {
 169         int i;
 170         int count;
 171         int temp;
 172         double val;
 173         int last_i;
 174
 175         last_i = 0;
 176
 177         if(win_size == 1)
 178         {
 179                 return(sorted_list[0]);
 180         }
 181
 182         /*
 183          * sweep through sorted list counting only those that are
 184          * within a win_size of now.  stop when we have win_size/2
 185          */
 186         count = 0;
 187         for(i=0; i < list_size; i++)
 188         {
 189                 if(sorted_ts[i] > (now - win_size))
 190                 {
 191                         count++;
 192                         /*
 193                          * need to remember this one for the even case
 194                          */
 195                         if(count == (win_size / 2))
 196                         {
 197                                 last_i = i;
 198                         }
 199                 }
 200
 201                 if(count == ((win_size / 2) + 1))
 202                         break;
 203         }
 204
 205         /*
 206          * i should now be the index of middle element -- almost
 207          */
 208         temp = win_size / 2;
 209         /*
 210          * if even
 211          */
 212         if(win_size == (temp * 2))
 213         {
 214                 val = (sorted_list[last_i] + sorted_list[i]) / 2.0;
 215         }
 216         else
 217         {
 218                 val = sorted_list[i];
 219         }
 220
 221         return(val);
 222 }
 223
 224 double
 225 FindTrimMedian(double *sorted_list,
 226            double *sorted_ts,
 227            double now,
 228            int win_size,
 229            double alpha)
 230 {
 231         int i;
 232         double val;
 233         int T;
 234         double vcount;
 235
 236         if(win_size == 1)
 237         {
 238                 return(sorted_list[0]);
 239         }
 240
 241         /*
 242          * sweep through sorted list counting only those that are
 243          * within a win_size of now.  stop when we have win_size/2
 244          */
 245         T = alpha * win_size;
 246         val = 0.0;
 247         vcount = 0.0;
 248         for(i=0; i < win_size ; i++)
 249         {
 250                 if((sorted_ts[i] > (now - win_size + T)) &&
 251                    (sorted_ts[i] < (now - T)))
 252                 {
 253                         val += sorted_list[i];
 254                         vcount++;
 255                 }
 256         }
 257
 258         if(vcount > 0.0)
 259         {
 260                 val = val / vcount;
 261         }
 262         else
 263         {
 264                 val = 0.0;
 265         }
 266
 267         return(val);
 268 }
 269
 270
 271 char *
 272 InitMedian(fbuff series, fbuff time_stamps, char *params)
 273 {
 274         struct median_state *state;
 275         char *p_str;
 276         double M_size;
 277         int i;
 278
 279         state = (struct median_state *)malloc(sizeof(struct median_state));
 280
 281         if(state == NULL)
 282         {
 283                 return(NULL);
 284         }
 285
 286         memset((void *)state,0,sizeof(struct median_state));
 287
 288         state->trim = 0;
 289         if(params != NULL)
 290         {
 291                 /*
 292                  * first parameter is siaze of median window
 293                  */
 294                 p_str = params;
 295                 M_size = (int)strtod(p_str,&p_str);
 296                 state->M_size = M_size;
 297         }
 298         else
 299         {
 300                 state->M_size = MAX_MED_SIZE;
 301         }
 302
 303         /*
 304          * init M_ts to all SMALL_TS so that MSort will build sorted
 305          * list correctly
 306          */
 307         for(i=0; i < MAX_MED_SIZE; i++)
 308         {
 309                 state->M_ts[i] = SMALL_TS;
 310         }
 311         state->artificial_time = 0;
 312         state->series = series;
 313         state->time_stamps = time_stamps;
 314         state->M_count = 0;
 315
 316         return((char *)state);
 317 }
 318
 319 void
 320 FreeMedian(char *state)
 321 {
 322         free(state);
 323         return;
 324 }
 325
 326 void
 327 UpdateMedian(char *state,
 328              double ts,
 329              double value)
 330 {
 331         struct median_state *s = (struct median_state *)state;
 332         int curr_size;
 333
 334
 335         curr_size = F_COUNT(s->series);
 336
 337         if(curr_size > s->M_size)
 338         {
 339                 s->M_count = curr_size = s->M_size;
 340         }
 341         else
 342         {
 343                 s->M_count = curr_size;
 344         }
 345
 346         /*
 347          * update the sorted list
 348          */
 349
 350         /*
 351          * increment the artificial time stamp
 352          */
 353         s->artificial_time = s->artificial_time + 1;
 354
 355         /*
 356          * use artificial time stamp instead of real one to
 357          * keep things in terms of entries instead of seconds
 358          */
 359         MSort(s->M_array,s->M_ts,value,s->artificial_time,curr_size);
 360
 361         return;
 362 }
 363
 364
 365 int
 366 ForcMedian(char *state, double *forecast)
 367 {
 368         struct median_state *s = (struct median_state *)state;
 369         int l_index;
 370         double val;
 371
 372         if(s->M_count < 1)
 373         {
 374                 return(0);
 375         }
 376
 377         l_index = s->M_count/2;
 378
 379         /*
 380          * if even
 381          */
 382         if(s->M_count == (2 * l_index))
 383         {
 384                 val = (s->M_array[l_index-1] + s->M_array[l_index]) / 2.0;
 385         }
 386         else
 387         {
 388                 val = s->M_array[l_index];
 389         }
 390
 391         *forecast = val;
 392         return(1);
 393 }
 394
 395 char *InitTrimMedian(fbuff series, fbuff time_stamps, char *params)
 396 {
 397         struct median_state *state;
 398         char *p_str;
 399         double M_size;
 400         double alpha;
 401         int i;
 402
 403         state = (struct median_state *)malloc(sizeof(struct median_state));
 404
 405         if(state == NULL)
 406         {
 407                 return(NULL);
 408         }
 409
 410         memset((void *)state,0,sizeof(struct median_state));
 411
 412         state->trim = 0;
 413         if(params != NULL)
 414         {
 415                 /*
 416                  * first parameter is size of median window
 417                  */
 418                 p_str = params;
 419                 M_size = (int)strtod(p_str,&p_str);
 420                 state->M_size = M_size;
 421
 422                 /*
 423                  * second parameter is the alpha trim value
 424                  */
 425                 alpha = strtod(p_str, &p_str);
 426                 state->trim = alpha;
 427         }
 428         else
 429         {
 430                 state->M_size = MAX_MED_SIZE;
 431                 state->trim = 0.0;
 432         }
 433
 434         /*
 435          * init M_ts to all SMALL_TS so that MSort will build sorted
 436          * list correctly
 437          */
 438         for(i=0; i < MAX_MED_SIZE; i++)
 439         {
 440                 state->M_ts[i] = SMALL_TS;
 441         }
 442         state->artificial_time = 0;
 443         state->series = series;
 444         state->time_stamps = time_stamps;
 445         state->M_count = 0;
 446
 447         return((char *)state);
 448 }
 449
 450 void
 451 FreeTrimMedian(char *state)
 452 {
 453         FreeMedian(state);
 454         return;
 455 }
 456
 457 void
 458 UpdateTrimMedian(char *state, double ts, double value)
 459 {
 460         UpdateMedian(state,ts,value);
 461         return;
 462 }
 463
 464
 465
 466 /*
 467  * page 264 in DSP book
 468  */
 469 int
 470 ForcTrimMedian(char *state, double *forecast)
 471 {
 472         struct median_state *s = (struct median_state *)state;
 473         double val;
 474         int l_index;
 475         int T;
 476
 477         if(s->M_count < 1)
 478         {
 479                 return(0);
 480         }
 481
 482
 483         T = (int)(s->trim * s->M_count);
 484
 485         val = 0.0;
 486         for(l_index = T; l_index < (s->M_count - T); l_index++)
 487         {
 488                 val += s->M_array[l_index];
 489         }
 490
 491         val = val / (double)(s->M_count - 2*T);
 492
 493
 494         *forecast = val;
 495
 496         return(1);
 497 }
 498
 499
 500
 501 char *InitAdMedian(fbuff series, fbuff time_stamps, char *params)
 502 {
 503         struct ad_median_state *state;
 504         char *p_str;
 505         double win;
 506         double min_win;
 507         double max_win;
 508         double offset;
 509         int i;
 510
 511         state = (struct ad_median_state *)
 512                 malloc(sizeof(struct ad_median_state));
 513
 514         if(state == NULL)
 515         {
 516                 return(NULL);
 517         }
 518
 519         memset((void *)state,0,sizeof(struct median_state));
 520
 521
 522         if(params != NULL)
 523         {
 524                 /*
 525                  * first parameter is initial target size of median window
 526                  */
 527                 p_str = params;
 528                 win = (int)strtod(p_str,&p_str);
 529                 state->win = win;
 530
 531                 /*
 532                  * defaults
 533                  */
 534                 min_win = win;
 535                 max_win = win;
 536                 offset = 1;
 537
 538                 /*
 539                  * second parameter is the min window value
 540                  */
 541                 min_win = strtod(p_str, &p_str);
 542
 543                 /*
 544                  * third parameter is the max window value
 545                  */
 546                 max_win = strtod(p_str, &p_str);
 547
 548                 /*
 549                  * fourth parameter is the offset value to use
 550                  * when calculating alternative window sizes
 551                  */
 552                 offset = strtod(p_str, &p_str);
 553                 state->offset = (int)offset;
 554
 555                 /*
 556                  * sanity checking
 557                  */
 558                 if(win == 1)
 559                         state->offset = 0;
 560
 561                 if((win - state->offset) < 0)
 562                         state->offset = 1;
 563
 564                 if(min_win > (win - state->offset))
 565                         min_win = win - state->offset;
 566
 567                 if(max_win < (win + state->offset))
 568                         max_win = win + state->offset;
 569
 570                 state->min = (int)min_win;
 571                 state->max = (int)max_win;
 572         }
 573         else
 574         {
 575                 state->min = state->max = state->win = MAX_MED_SIZE;
 576                 state->offset = 1;
 577         }
 578
 579         /*
 580          * init M_ts to all SMALL_TS so that MSort will build sorted
 581          * list correctly
 582          */
 583         for(i=0; i < MAX_MED_SIZE; i++)
 584         {
 585                 state->M_ts[i] = SMALL_TS;
 586         }
 587         state->artificial_time = 0;
 588         state->series = series;
 589         state->time_stamps = time_stamps;
 590         state->M_count = 0;
 591         state->artificial_time = 1;
 592
 593         return((char *)state);
 594 }
 595
 596 void
 597 FreeAdMedian(char *state)
 598 {
 599                 free(state);
 600                 return;
 601 }
 602
 603 void
 604 UpdateAdMedian(char *state, double ts, double value)
 605 {
 606         struct ad_median_state *s = (struct ad_median_state *)state;
 607         int curr_size;
 608         int win;
 609         double less_val;
 610         double eq_val;
 611         double more_val;
 612         double less_err;
 613         double eq_err;
 614         double more_err;
 615         int lo_offset;
 616         int hi_offset;
 617
 618
 619         curr_size = F_COUNT(s->series);
 620
 621         /*
 622          * M_size is the current window size, and M_count is the
 623          * current amount of data in the median buffer
 624          */
 625
 626         if(curr_size > s->max)
 627         {
 628                 s->M_count = curr_size = s->max;
 629         }
 630         else
 631         {
 632                 s->M_count = curr_size;
 633         }
 634
 635         /*
 636          * update the sorted list
 637          */
 638
 639         /*
 640          * increment the artificial time stamp
 641          */
 642         s->artificial_time = s->artificial_time + 1;
 643
 644         /*
 645          * use artificial time stamp instead of real one to
 646          * keep things in terms of entries instead of seconds
 647          */
 648         MSort(s->M_array,s->M_ts,value,s->artificial_time,curr_size);
 649
 650
 651         /*
 652          * calculate the window based on how much data there is
 653          */
 654         if(curr_size > s->win)
 655         {
 656                 win = s->win;
 657         }
 658         else
 659         {
 660                 win = curr_size;
 661         }
 662         /*
 663          * find the median using the current
 664          * window size
 665          */
 666         eq_val = FindMedian(s->M_array,
 667                             s->M_ts,
 668                             s->M_count,
 669                             s->artificial_time,
 670                             win);
 671
 672         /*
 673          * we want to wait until there is enough data before we start
 674          * to adapt.  We don't start to adjust s->win until there is
 675          * enough data to get out to the max window size
 676          */
 677         if(curr_size < s->max)
 678         {
 679                 return;
 680         }
 681
 682         if((win - s->offset) < 0)
 683         {
 684                 lo_offset = win - 1;
 685         }
 686         else
 687         {
 688                 lo_offset = s->offset;
 689         }
 690
 691         if((win + s->offset) > s->M_count)
 692         {
 693                 hi_offset = s->M_count - win - 1;
 694         }
 695         else
 696         {
 697                 hi_offset = s->offset;
 698         }
 699
 700         /*
 701          * find the median for a smaller window -- offset
 702          * controls how much smaller or bigger the window should be
 703          * that we consider
 704          */
 705         less_val = FindMedian(s->M_array,
 706                               s->M_ts,
 707                               s->M_count,
 708                               s->artificial_time,
 709                               lo_offset);
 710
 711         /*
 712          * find the median for a bigger window -- offset
 713          * controls how much smaller or bigger the window should be
 714          * that we consider
 715          */
 716         more_val = FindMedian(s->M_array,
 717                               s->M_ts,
 718                               s->M_count,
 719                               s->artificial_time,
 720                               hi_offset);
 721
 722         /*
 723          * now, calculate the errors
 724          */
 725         less_err = (value - less_val) * (value - less_val);
 726         more_err = (value - more_val) * (value - more_val);
 727         eq_err = (value - eq_val) * (value - eq_val);
 728
 729         /*
 730          * adapt the window according to the direction giving us the
 731          * smallest error
 732          */
 733         if(less_err < eq_err)
 734         {
 735                 if(less_err < more_err)
 736                 {
 737                         win = win - 1;
 738                 }
 739                 else if(more_err < eq_err)
 740                 {
 741                         win = win + 1;
 742                 }
 743         }
 744         else if(more_err < eq_err)
 745         {
 746                 if(more_err < less_err)
 747                 {
 748                         win = win + 1;
 749                 }
 750                 else if(less_err < eq_err)
 751                 {
 752                         win = win - 1;
 753                 }
 754         }
 755
 756         s->win = win;
 757
 758         return;
 759 }
 760
 761
 762 int
 763 ForcAdMedian(char *state, double *forecast)
 764 {
 765         struct ad_median_state *s = (struct ad_median_state *)state;
 766         double forc;
 767
 768         if(s->M_count == 0)
 769         {
 770                 return(0);
 771         }
 772
 773         if(s->M_count < s->max)
 774         {
 775                 forc = FindMedian(s->M_array,
 776                                   s->M_ts,
 777                                   s->M_count,
 778                                   s->artificial_time,
 779                                   s->M_count);
 780         }
 781         else
 782         {
 783                 forc = FindMedian(s->M_array,
 784                                   s->M_ts,
 785                                   s->M_count,
 786                                   s->artificial_time,
 787                                   s->win);
 788
 789         }
 790
 791         *forecast = forc;
 792
 793         return(1);
 794 }
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804 #ifdef NOTYET
 805
 806 AdjustedMedian(tp,pred,min,max)
 807 struct thruput_pred *tp;
 808 int pred;
 809 int min;
 810 int max;
 811 {
 812         double err1;
 813         double err2;
 814         double err3;
 815         double val1;
 816         double val2;
 817         double val3;
 818         int win;
 819         struct value_el *data = tp->data;
 820         int prev = MODMINUS(tp->head,1,HISTORYSIZE);
 821         int head = tp->head;
 822         double value = tp->data[head].value;
 823         double predictor = PRED(pred,data[prev]);
 824         int out_win;
 825         double out_val;
 826
 827         win = MedWIN(pred,data[prev]);
 828         if(MODMINUS(head,tp->tail,HISTORYSIZE) < 3)
 829         {
 830                 PRED(pred,data[head]) = value;
 831                 MedWIN(pred,data[head]) = win;
 832                 return(0.0);
 833         }
 834         if(win < 2)
 835                 win = 2;
 836         val1 = Median(tp,win-1);
 837         val2 = Median(tp,win+1);
 838         val3 = Median(tp,win);
 839
 840         err1 = (value - val1) * (value - val1);
 841         err2 = (value - val2) * (value - val2);
 842         err3 = (value - val3) * (value - val3);
 843
 844         if(err1 < err2)
 845         {
 846                 if(err1 < err3)
 847                 {
 848                         out_win = win - 1;
 849                         out_val = val1;
 850                 }
 851                 else
 852                 {
 853                         out_win = win;
 854                         out_val = val3;
 855                 }
 856         }
 857         else if(err2 < err1)
 858         {
 859                 if(err2 < err3)
 860                 {
 861                         out_win = win + 1;
 862                         out_val = val2;
 863                 }
 864                 else
 865                 {
 866                         out_win = win;
 867                         out_val = val3;
 868                 }
 869         }
 870         else
 871         {
 872                 out_win = win;
 873                 out_val = val3;
 874         }
 875
 876         if(out_win < min)
 877         {
 878                 out_win = min;
 879                 out_val = Median(tp,min);
 880         }
 881
 882         if(out_win > max)
 883         {
 884                 out_win = max;
 885                 out_val = Median(tp,max);
 886         }
 887
 888         MedWIN(pred,data[head]) = out_win;
 889         PRED(pred,data[head]) = out_val;
 890 }
 891
 892 AdjustedMean(tp,pred,min,max)
 893 struct thruput_pred *tp;
 894 int pred;
 895 int min;
 896 int max;
 897 {
 898         double err1;
 899         double err2;
 900         double err3;
 901         double val1;
 902         double val2;
 903         double val3;
 904         int win;
 905         struct value_el *data = tp->data;
 906         int prev = MODMINUS(tp->head,1,HISTORYSIZE);
 907         int head = tp->head;
 908         double value = tp->data[head].value;
 909         double predictor = PRED(pred,data[prev]);
 910         int out_win;
 911         double out_val;
 912
 913         win = MeanWIN(pred,data[prev]);
 914         if(MODMINUS(head,tp->tail,HISTORYSIZE) < 3)
 915         {
 916                 PRED(pred,data[head]) = value;
 917                 MeanWIN(pred,data[head]) = WIN_INIT;
 918                 return(0.0);
 919         }
 920         if(win < 2)
 921                 win = 2;
 922         val1 = TrimmedMedian(tp,win-1,0.0);
 923         val2 = TrimmedMedian(tp,win+1,0.0);
 924         val3 = TrimmedMedian(tp,win,0.0);
 925
 926         err1 = (value - val1) * (value - val1);
 927         err2 = (value - val2) * (value - val2);
 928         err3 = (value - val3) * (value - val3);
 929
 930         if(err1 < err2)
 931         {
 932                 if(err1 < err3)
 933                 {
 934                         out_win = win - 1;
 935                         out_val = val1;
 936                 }
 937                 else
 938                 {
 939                         out_win = win;
 940                         out_val = val3;
 941                 }
 942         }
 943         else if(err2 < err1)
 944         {
 945                 if(err2 < err3)
 946                 {
 947                         out_win = win + 1;
 948                         out_val = val2;
 949                 }
 950                 else
 951                 {
 952                         out_win = win;
 953                         out_val = val3;
 954                 }
 955         }
 956         else
 957         {
 958                 out_win = win;
 959                 out_val = val3;
 960         }
 961
 962         if(out_win < min)
 963         {
 964                 out_win = min;
 965                 out_val = TrimmedMedian(tp,min,0.0);
 966         }
 967
 968         if(out_win > max)
 969         {
 970                 out_win = max;
 971                 out_val = TrimmedMedian(tp,max,0.0);
 972         }
 973
 974         MeanWIN(pred,data[head]) = out_win;
 975         PRED(pred,data[head]) = out_val;
 976 }
 977
 978 #endif