00001 #ifndef PARTICLEFILTER_H
00002 #define PARTICLEFILTER_H
00003 #include <stdlib.h>
00004 #include<vector>
00005 #include<utility>
00006 #include<utils/gvalues.h>
00007
00008
00009
00020 typedef std::pair<uint,uint> UIntPair;
00021
00022 template <class OutputIterator, class Iterator>
00023 double toNormalForm(OutputIterator& out, const Iterator & begin, const Iterator & end){
00024
00025 double lmax=-MAXDOUBLE;
00026 for (Iterator it=begin; it!=end; it++){
00027 lmax=lmax>((double)(*it))? lmax: (double)(*it);
00028 }
00029
00030 for (Iterator it=begin; it!=end; it++){
00031 *out=exp((double)(*it)-lmax);
00032 out++;
00033 }
00034 return lmax;
00035 }
00036
00037 template <class OutputIterator, class Iterator, class Numeric>
00038 void toLogForm(OutputIterator& out, const Iterator & begin, const Iterator & end, Numeric lmax){
00039
00040 for (Iterator it=begin; it!=end; it++){
00041 *out=log((Numeric)(*it))-lmax;
00042 out++;
00043 }
00044 return lmax;
00045 }
00046
00047 template <class WeightVector>
00048 void resample(std::vector<int>& indexes, const WeightVector& weights, unsigned int nparticles=0){
00049 double cweight=0;
00050
00051
00052 unsigned int n=0;
00053 for (typename WeightVector::const_iterator it=weights.begin(); it!=weights.end(); ++it){
00054 cweight+=(double)*it;
00055 n++;
00056 }
00057
00058 if (nparticles>0)
00059 n=nparticles;
00060
00061
00062 double interval=cweight/n;
00063
00064
00065 double target=interval*::drand48();
00066
00067
00068 cweight=0;
00069 indexes.resize(n);
00070
00071 n=0;
00072 unsigned int i=0;
00073 for (typename WeightVector::const_iterator it=weights.begin(); it!=weights.end(); ++it, ++i){
00074 cweight+=(double)* it;
00075 while(cweight>target){
00076 indexes[n++]=i;
00077 target+=interval;
00078 }
00079 }
00080 }
00081
00082 template <typename WeightVector>
00083 void normalizeWeights(WeightVector& weights, unsigned int size, double minWeight){
00084 double wmin=MAXDOUBLE;
00085 double wmax=-MAXDOUBLE;
00086 for (uint i=0; i<size; i++){
00087 wmin=wmin<weights[i]?wmin:weights[i];
00088 wmax=wmax>weights[i]?wmax:weights[i];
00089 }
00090 double min_normalized_value=log(minWeight);
00091 double max_normalized_value=log(1.);
00092 double dn=max_normalized_value-min_normalized_value;
00093 double dw=wmax-wmin;
00094 if (dw==0) dw=1;
00095 double scale=dn/dw;
00096 double offset=-wmax*scale;
00097 for (uint i=0; i<size; i++){
00098 double w=weights[i];
00099 w=scale*w+offset;
00100 weights[i]=exp(w);
00101 }
00102 }
00103
00104 template <typename Vector>
00105 void repeatIndexes(Vector& dest, const std::vector<int>& indexes, const Vector& particles){
00106
00107
00108
00109
00110
00111
00112
00113 dest.resize(indexes.size());
00114
00115 unsigned int i=0;
00116 for (std::vector<int>::const_iterator it=indexes.begin(); it!=indexes.end(); ++it){
00117 dest[i]=particles[*it];
00118 i++;
00119 }
00120 }
00121
00122 template <typename Vector>
00123 void repeatIndexes(Vector& dest, const std::vector<int>& indexes2, const Vector& particles, const std::vector<int>& indexes){
00124
00125 dest=particles;
00126 unsigned int i=0;
00127 for (std::vector<int>::const_iterator it=indexes2.begin(); it!=indexes2.end(); ++it){
00128 dest[indexes[i]]=particles[*it];
00129 i++;
00130 }
00131 }
00132
00133
00134 template <class Iterator>
00135 double neff(const Iterator& begin, const Iterator& end){
00136 double sum=0;
00137 for (Iterator it=begin; it!=end; ++it){
00138 sum+=*it;
00139 }
00140 double cum=0;
00141 for (Iterator it=begin; it!=end; ++it){
00142 double w=*it/sum;
00143 cum+=w*w;
00144 }
00145 return 1./cum;
00146 }
00147
00148
00149
00150 template <class OutputIterator, class Iterator>
00151 void rle(OutputIterator& out, const Iterator & begin, const Iterator & end){
00152 unsigned int current=0;
00153 unsigned int count=0;
00154 for (Iterator it=begin; it!=end; it++){
00155 if (it==begin){
00156 current=*it;
00157 count=1;
00158 continue;
00159 }
00160 if (((uint)*it) ==current)
00161 count++;
00162 if (((uint)*it)!=current){
00163 *out=std::make_pair(current,count);
00164 out++;
00165 current=*it;
00166 count=1;
00167 }
00168 }
00169 if (count>0)
00170 *out=std::make_pair(current,count);
00171 out++;
00172 }
00173
00174 #endif
00175