00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 
00014 
00015 
00016 
00017 
00018 
00019 
00020 
00021 
00022 
00023 
00024 
00025 
00026 
00027 
00028 
00029 
00030 
00031 
00032 
00033 
00034 
00035 #ifndef CEVENT_COUNTER_DATA_H
00036 #define CEVENT_COUNTER_DATA_H
00037 #include<cmath>
00038 #include <iostream>
00039 #include <cstdio>
00040 #include <Eigen/Dense>
00041 #include <inttypes.h>
00042 
00044 #define EVENTMAP_OCCU  255
00045 #define EVENTMAP_UNKNOWN  127
00046 #define EVENTMAP_FREE  0
00047 #define EVENTMAP_NUMBER_OF_EVENTS_STORED 128
00048 
00049 
00051 #define EVENTMAP_USE_RECENCY_FILTERING
00052 #define EVENTMAP_OBSERVATION_LIMIT 30000.0
00053 #define Wforget (5000.0/5001.0)
00054 
00058 struct TEventData
00059 {
00060     uint8_t     occval;                         
00061     float               a_exit_event; 
00062     float               b_exit_event; 
00063 
00064     float               a_entry_event;
00065     float               b_entry_event; 
00066     uint64_t  events;        
00067 
00068     TEventData():
00069         occval(127),
00070         a_exit_event(1.0),
00071         b_exit_event(1.0),
00072         a_entry_event(1.0),
00073         b_entry_event(1.0),
00074         events(0)
00075     {
00076     };
00077 
00078     ~TEventData()
00079     {
00080     }
00081 
00082 
00083     const TEventData &operator =(const TEventData ©)
00084     {
00085         occval = copy.occval;
00086         a_exit_event = copy.a_exit_event;
00087         b_exit_event = copy.b_exit_event;
00088         a_entry_event = copy.a_entry_event;
00089         b_entry_event = copy.b_entry_event;
00090         events = copy.events;
00091 
00092         return *this;
00093     }
00094 
00095     TEventData(const TEventData& copy):
00096         occval(copy.occval),
00097         a_exit_event(copy.a_exit_event),
00098         b_exit_event(copy.b_exit_event),
00099         a_entry_event(copy.a_entry_event),
00100         b_entry_event(copy.b_entry_event),
00101         events( copy.events)
00102     {
00103     }
00104 
00105 
00106     void Reset()
00107     {
00108         occval = 127;
00109         a_exit_event = 1;
00110         b_exit_event = 1;
00111         a_entry_event = 1;
00112         b_entry_event = 1;
00113     }
00114 
00115     bool getBit(int ind)
00116     {
00117         if((events & (1<<ind))>0) return true;
00118         return false;
00119     }
00120 
00121 
00122     float computeShortTermOccupancy()
00123     {
00124         static float OCC_LOG_PROP = log(0.6/0.4);
00125         static float EMP_LOG_PROP = log(0.3/0.7);
00126         uint bs=0;
00127 
00128         if( (b_entry_event+b_exit_event-2) > 63) bs=64;
00129         else bs = b_entry_event+b_exit_event-2;
00130         uint cnt=0;
00131         float log_ogg_prob  = 0 ;
00132         for(uint i = 0; i<bs; i++)
00133         {
00134             if(getBit(i))
00135             {
00136                 log_ogg_prob += OCC_LOG_PROP;
00137                 cnt++;
00138             }
00139             else
00140             {
00141                 log_ogg_prob += EMP_LOG_PROP;
00142             }
00143         }
00144         return (1.0- 1.0/(1.0+exp(log_ogg_prob)));
00145     }
00146 
00147     int getObservations()
00148     {
00149         return (b_entry_event+b_exit_event);
00150     }
00151 
00155     void updateSimple(uint8_t meas)
00156     {
00157         events = events<<1;
00158         if(meas > EVENTMAP_UNKNOWN)  
00159         {
00160             events |= 1;
00161         }
00162         if(occval > EVENTMAP_UNKNOWN)  
00163         {
00164             updateExitEvent(meas, 1.0);
00165         }
00166         else if(occval < EVENTMAP_UNKNOWN)
00167         {
00168             updateEntryEvent(meas,1.0);
00169         }
00170         else
00171         {
00172             occval = meas; 
00173         }
00174     }
00175 
00179     void updateSimple(uint8_t meas, float w)
00180     {
00181         events = events<<1;
00182         if(meas > EVENTMAP_UNKNOWN)  
00183         {
00184             events |= 1;
00185         }
00186         if(occval > EVENTMAP_UNKNOWN)  
00187         {
00188             updateExitEvent(meas, w);
00189         }
00190         else if(occval < EVENTMAP_UNKNOWN)
00191         {
00192             updateEntryEvent(meas,w);
00193         }
00194         else
00195         {
00196             occval = meas; 
00197         }
00198 #ifdef EVENTMAP_USE_RECENCY_FILTERING
00199         performRecencyFiltering();
00200 #endif
00201     }
00202 
00203     void performRecencyFiltering()
00204     {
00205         if(occval > EVENTMAP_UNKNOWN)  
00206         {
00207             if(b_exit_event>EVENTMAP_OBSERVATION_LIMIT)
00208             {
00209                 float w = EVENTMAP_OBSERVATION_LIMIT / b_exit_event;
00210                 b_exit_event *= w;
00211                 a_exit_event *= w;
00212             }
00213             a_entry_event = 1.0 + (a_entry_event-1.0) * Wforget;
00214             b_entry_event = 1.0 + (b_entry_event-1.0) * Wforget;
00215         }
00216         else
00217         {
00218             if(b_entry_event>EVENTMAP_OBSERVATION_LIMIT)
00219             {
00220                 float w = EVENTMAP_OBSERVATION_LIMIT / b_entry_event;
00221                 b_entry_event *= w;
00222                 a_entry_event *= w;
00223             }
00224             a_exit_event = 1.0 + (a_exit_event-1.0) * Wforget;
00225             b_exit_event = 1.0 + (b_exit_event-1.0) * Wforget;
00226         }
00227     }
00228 
00229 
00233     void updateExitEvent(uint8_t meas, float w)
00234     {
00235         if(meas>EVENTMAP_UNKNOWN)  
00236         {
00237             b_exit_event+=1.0;
00238         }
00239         else
00240         {
00241             if(w>1.0) fprintf(stderr,"EXIT W larger that one = %f\n ",w);
00242             a_exit_event+=w;
00243             b_exit_event+=1.0;
00244             occval = EVENTMAP_FREE;
00245         }
00246     }
00250     void updateEntryEvent(uint8_t meas, float w)
00251     {
00252         if(meas>EVENTMAP_UNKNOWN)  
00253         {
00254             b_entry_event+=1.0;
00255             if(w>1.0) fprintf(stderr,"W larger that one = %f\n ",w);
00256             a_entry_event+=w;
00257             occval = EVENTMAP_OCCU;
00258         }
00259         else
00260         {
00261             b_entry_event+=1.0;
00262         }
00263     }
00264 
00266     int getEntryN()
00267     {
00268         return (int)((float)b_entry_event/(float)a_entry_event +0.5);
00269     }
00271     int getExitN()
00272     {
00273         return (int)((float)b_exit_event/(float)a_exit_event +0.5);
00274     }
00275 
00276 
00277     double entryL()
00278     {
00279         if(b_entry_event == 0) fprintf(stderr,"B Zero, which is not possible!!\n");
00280         return ((double)a_entry_event/(double)b_entry_event );
00281     }
00282     double exitL()
00283     {
00284         if(b_exit_event == 0) fprintf(stderr,"B Zero, which is not possible!!\n");
00285         return ((double)a_exit_event/(double)b_exit_event);
00286     }
00287     double L()
00288     {
00289         if(b_exit_event == 0) fprintf(stderr,"B Zero, which is not possible!!\n");
00290         return ((double)(a_exit_event+a_entry_event)/(double)(b_exit_event+b_entry_event));
00291     }
00292 
00293     double fac(int n)
00294     {
00295         double t=1;
00296         for (int i=n; i>1; i--)
00297             t*=i;
00298         return t;
00299     }
00300 
00301     double Bin(int n,double p,int r)
00302     {
00303         return fac(n)/(fac(n-r)*fac(r))*pow(p,r)*pow(1-p,n-r);
00304     }
00305 
00306     float binaryBayesUpdate(float Pold, float P)
00307     {
00308         return (( Pold*P) / ( Pold*P + (1.0f-Pold)*(1.0f-P)));  
00309     }
00310 
00314     void normalizeProb(float &p)
00315     {
00316         if(p>0.999) p = 0.999;
00317         if(p<0.001) p = 0.001;
00318 
00319     }
00324 
00327     float getOccStaticLikelihood()
00328     {
00329         if( (b_entry_event + b_exit_event)<20.0 )
00330         {
00331             return 0.5f;
00332         }
00333         Eigen::Matrix2f P;
00334         Eigen::Vector2f u1(0.5, 0.5);
00335         Eigen::Vector2f P0;
00336         float Lex = exitL();
00337         float Len = entryL();
00338         P(0,0) = (1.0-Len);
00339         P(0,1) =  Len;
00340         P(1,0)   = Lex;
00341         P(1,1) = (1-Lex);
00342 
00343         P0 = u1.transpose() * P;
00344         return (P0(1));
00345     }
00346 
00350     float getFreeStaticLikelihood()
00351     {
00352         if( (b_entry_event + b_exit_event)<20.0 )
00353         {
00354             return 0.5f;
00355         }
00356         Eigen::Matrix2f P;
00357         Eigen::Vector2f P0;
00358         Eigen::Vector2f u1(0.5, 0.5);
00359         float Lex = exitL();
00360         float Len = entryL();
00361         P(0,0) = (1.0-Len);
00362         P(0,1) =  Len;
00363         P(1,0)   = Lex;
00364         P(1,1) = (1-Lex);
00365 
00366         P0 = u1.transpose() * P;
00367         return (P0(0));
00368     }
00369 
00377     float computeSemiStaticLikelihood(int N)
00378     {
00379         if( (b_entry_event + b_exit_event)<20.0 )
00380         {
00381             return 0.0f;
00382         }
00383 
00384 
00385         Eigen::Matrix2f P;
00386         Eigen::Vector2f u2(0, 1.0);
00387         Eigen::Vector2f u3(1.0, 0);
00388         float Lex = exitL();
00389         float Len = entryL();
00390         Eigen::Vector2f P2,P3;
00391         P(0,0) = (1.0-Len);
00392         P(0,1) =  Len;
00393         P(1,0)   = Lex;
00394         P(1,1) = (1-Lex);
00395         for(int i=0; i<N; i++) P = P*P;
00396         P2 = u2.transpose() * P;
00397         P3 = u3.transpose() * P;
00398 
00399         float Po = P2(1);
00400         float Pu = binaryBayesUpdate(Po, P3(0));
00401         normalizeProb(Pu);
00402 
00403         return Pu;
00404     }
00405 
00406     float getOccupancyNow()
00407     {
00408         if( (b_entry_event + b_exit_event)<50.0)
00409         {
00410             return 0.5f;
00411         }
00412         float Po = 0;
00413 
00414 
00415         Po = computeShortTermOccupancy();
00416         float Lex = exitL();
00417         float Len = entryL();
00418         Eigen::Matrix2f P;
00419         P(0,0) = (1.0-Len);
00420         P(0,1) =  Len;
00421         P(1,0)   = Lex;
00422         P(1,1) = (1-Lex);
00423         Eigen::Vector2f u1(1.0-Po, Po);
00424         Eigen::Vector2f u = u1.transpose() *P;
00425         return u(1);
00426     }
00427 
00433 
00434     float predictOccupancy(int N)
00435     {
00436         if( (b_entry_event + b_exit_event)<20.0)
00437         {
00438             return 0.5f;
00439         }
00441         float ps = getOccStaticLikelihood();
00442         if(ps>0.6) return 0.9;
00444 
00445 
00446         float Po = computeShortTermOccupancy();
00447         float Lex = exitL();
00448         float Len = entryL();
00449         Eigen::Matrix2f P;
00450         P(0,0) = (1.0-Len);
00451         P(0,1) =  Len;
00452         P(1,0)   = Lex;
00453         P(1,1) = (1-Lex);
00454         Eigen::Vector2f u1(1.0-Po, Po);
00455         for(int i=0; i<N; i++) P = P*P;
00456         Eigen::Vector2f u = u1.transpose() *P;
00457 
00458         return u(1);
00459     }
00460 
00461 
00462 };
00463 #endif
00464