discretepdf.cpp
Go to the documentation of this file.
00001 // $Id: discretepdf.cpp 29890 2009-02-02 10:22:01Z tdelaet $
00002 // Copyright (C) 2002 Klaas Gadeyne <first dot last at gmail dot com>
00003 //
00004 // This program is free software; you can redistribute it and/or modify
00005 // it under the terms of the GNU Lesser General Public License as published by
00006 // the Free Software Foundation; either version 2.1 of the License, or
00007 // (at your option) any later version.
00008 //
00009 // This program is distributed in the hope that it will be useful,
00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 // GNU Lesser General Public License for more details.
00013 //
00014 // You should have received a copy of the GNU Lesser General Public License
00015 // along with this program; if not, write to the Free Software
00016 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00017 //
00018 #include "discretepdf.h"
00019 #include "../bfl_err.h"
00020 #include "../wrappers/rng/rng.h"
00021 #include <vector>
00022 #include <iostream>
00023 
00024 
00025 namespace BFL
00026 {
00027   using namespace std;
00028   using namespace MatrixWrapper;
00029 
00030 
00031   DiscretePdf::DiscretePdf(unsigned int num_states): Pdf<int>(1)
00032         ,_num_states(num_states)
00033   {
00034     //discrete pdf has dimension 1
00035     _Values_p = new vector<Probability>(num_states);
00036     for (int i=0; i<NumStatesGet();i++)
00037     {
00038         (*_Values_p)[i] = (Probability)(1.0/NumStatesGet());
00039     }
00040     _CumPDF.insert(_CumPDF.begin(),num_states+1,0.0);
00041     CumPDFUpdate();
00042 #ifdef __CONSTRUCTOR__
00043     cout << "DiscretePdf constructor\n";
00044 #endif // __CONSTRUCTOR__
00045   }
00046 
00047   DiscretePdf::DiscretePdf(const DiscretePdf & my_dpdf):Pdf<int>(my_dpdf)
00048         ,_num_states(my_dpdf.NumStatesGet())
00049   {
00050     _Values_p = new vector<Probability>(this->NumStatesGet());
00051     (*_Values_p) = my_dpdf.ProbabilitiesGet();
00052     _CumPDF.insert(_CumPDF.begin(),NumStatesGet()+1,0.0);
00053     CumPDFUpdate();
00054 #ifdef __CONSTRUCTOR__
00055     cout << "DiscretePdf copy constructor\n";
00056 #endif // __CONSTRUCTOR__
00057   }
00058 
00059   DiscretePdf::~DiscretePdf()
00060   {
00061 #ifdef __CONSTRUCTOR__
00062     cout << "DiscretePdf destructor\n";
00063 #endif
00064     // Release memory!
00065     delete _Values_p;
00066   }
00067 
00068   //Clone function
00069   DiscretePdf* DiscretePdf::Clone() const
00070   {
00071       return new DiscretePdf(*this);
00072   }
00073 
00074   unsigned int DiscretePdf::NumStatesGet()const
00075   {
00076     return _num_states;
00077   }
00078 
00079 
00080   Probability DiscretePdf::ProbabilityGet(const int& state) const
00081   {
00082     assert((int)state >= 0 && state < NumStatesGet());
00083     return (*_Values_p)[state];
00084   }
00085 
00086   bool DiscretePdf::ProbabilitySet(int state, Probability a)
00087   {
00088     assert((int)state >= 0 && state < NumStatesGet());
00089     assert(a<=1);
00090 
00091     // renormalize other probabilities such that sum of probabilities will be
00092     // one after the probability of state is set to a
00093     // This should keep the probabilities normalized
00094     Probability old_prob_state = ProbabilityGet(state);
00095     if (old_prob_state!=1) {
00096         double normalization_factor = (1-a)/(1-old_prob_state);
00097         for (int i=0; i<NumStatesGet();i++)
00098         {
00099             (*_Values_p)[i] = (Probability)( (double)( (*_Values_p)[i] )* normalization_factor);
00100         }
00101     }
00102     else{
00103         for (int i=0; i<NumStatesGet();i++)
00104         {
00105             (*_Values_p)[i] = (Probability)( (1-a)/(NumStatesGet()-1));
00106         }
00107     }
00108     (*_Values_p)[state] = a;
00109     return CumPDFUpdate();
00110   }
00111 
00112   vector<Probability> DiscretePdf::ProbabilitiesGet() const
00113   {
00114     return *_Values_p;
00115   }
00116 
00117   bool DiscretePdf::ProbabilitiesSet(vector<Probability> & v)
00118   {
00119     assert(v.size() == NumStatesGet());
00120 
00121     *_Values_p = v;
00122     //normalize the probabilities and update the cumulative sum
00123     return (NormalizeProbs() && CumPDFUpdate());
00124   }
00125 
00126   // For optimal performance!
00127   bool
00128   DiscretePdf::SampleFrom (vector<Sample<int> >& list_samples,
00129                            const unsigned int num_samples,
00130                            int method,
00131                            void * args) const
00132   {
00133     switch(method)
00134     {
00135       case DEFAULT: // O(N log(N) efficiency)
00136           return Pdf<int>::SampleFrom(list_samples, num_samples,method,args);
00137 
00138       case RIPLEY: // See mcpdf.cpp for more explanation
00139           {
00140             list_samples.resize(num_samples);
00141             // GENERATE N ORDERED IID UNIFORM SAMPLES
00142             std::vector<double> unif_samples(num_samples);
00143             for ( unsigned int i = 0; i < num_samples ; i++)
00144               unif_samples[i] = runif();
00145 
00146             /* take n-th racine of u_N */
00147             unif_samples[num_samples-1] = pow(unif_samples[num_samples-1],
00148                                            double (1.0/num_samples));
00149             /* rescale samples */
00150             for ( int i = num_samples-2; i >= 0 ; i--)
00151                 unif_samples[i] = pow(unif_samples[i], double (1.0/(i+1))) * unif_samples[i+1];
00152 
00153             // CHECK WHERE THESE SAMPLES ARE IN _CUMPDF
00154             unsigned int index = 0;
00155             unsigned int num_states = NumStatesGet();
00156             vector<double>::const_iterator CumPDFit = _CumPDF.begin();
00157             vector<Sample<int> >::iterator sit = list_samples.begin();
00158 
00159             for ( unsigned int i = 0; i < num_samples ; i++)
00160               {
00161                 while ( unif_samples[i] > *CumPDFit )
00162                 {
00163                 // check for internal error
00164                 assert(index <= num_states);
00165                 index++; CumPDFit++;
00166                 }
00167           int a = index - 1;
00168               sit->ValueSet(a);
00169               sit++;
00170               }
00171             return true;
00172           }
00173       default:
00174             cerr << "DiscretePdf::Samplefrom(int, void *): No such sampling method" << endl;
00175             return false;
00176     }
00177   }
00178 
00179 
00180 
00181   bool DiscretePdf::SampleFrom (Sample<int>& one_sample, int method, void * args) const
00182   {
00183     switch(method)
00184       {
00185       case DEFAULT:
00186         {
00187           // Sample from univariate uniform rng between 0 and 1;
00188           double unif_sample; unif_sample = runif();
00189           // Compare where we should be
00190           unsigned int index = 0;
00191           while ( unif_sample > _CumPDF[index] )
00192             {
00193               assert(index <= NumStatesGet());
00194               index++;
00195             }
00196           int a = index - 1;
00197           one_sample.ValueSet(a);
00198           return true;
00199         }
00200       default:
00201         cerr << "DiscretePdf::Samplefrom(int, void *): No such sampling method"
00202              << endl;
00203         return false;
00204       }
00205   }
00206 
00207   bool DiscretePdf::NormalizeProbs()
00208   {
00209     double SumOfProbs = 0.0;
00210     for ( unsigned int i = 0; i < NumStatesGet() ; i++){
00211         SumOfProbs += (*_Values_p)[i];
00212     }
00213     if (SumOfProbs > 0){
00214       for ( unsigned int i = 0; i < NumStatesGet() ; i++){
00215           (*_Values_p)[i] = (Probability)( (double) ( (*_Values_p)[i]) /SumOfProbs);
00216       }
00217       return true;
00218     }
00219     else{
00220       cerr << "DiscretePdf::NormalizeProbs(): SumOfProbs = " << SumOfProbs << endl;
00221       return false;
00222     }
00223   }
00224 
00225   bool DiscretePdf::CumPDFUpdate()
00226   {
00227     // precondition: sum of probabilities should be 1
00228     double CumSum=0.0;
00229     static vector<double>::iterator CumPDFit;
00230     CumPDFit = _CumPDF.begin();
00231     *CumPDFit = 0.0;
00232 
00233     // Calculate the Cumulative PDF
00234     for ( unsigned int i = 0; i < NumStatesGet(); i++)
00235     {
00236             CumPDFit++;
00237             // Calculate the __normalised__ Cumulative sum!!!
00238             CumSum += ( (*_Values_p)[i] );
00239             *CumPDFit = CumSum;
00240     }
00241     // Check if last element of valuelist is +- 1
00242     assert( (_CumPDF[NumStatesGet()] >= 1.0 - NUMERIC_PRECISION) &&
00243             (_CumPDF[NumStatesGet()] <= 1.0 + NUMERIC_PRECISION) );
00244 
00245     _CumPDF[NumStatesGet()]=1;
00246 
00247     return true;
00248   }
00249 
00250   int DiscretePdf::MostProbableStateGet()
00251   {
00252     int index_mostProbableState = -1;
00253     Probability prob_mostProbableState = 0.0;
00254     for (int state = 0 ; state < NumStatesGet() ; state++)
00255     {
00256        if ( (*_Values_p)[state] >= prob_mostProbableState)
00257        {
00258             index_mostProbableState = state;
00259             prob_mostProbableState = (*_Values_p)[state];
00260        }
00261     }
00262     return index_mostProbableState;
00263   }
00264 
00265 
00266 } // End namespace


bfl
Author(s): Klaas Gadeyne, Wim Meeussen, Tinne Delaet and many others. See web page for a full contributor list. ROS package maintained by Wim Meeussen.
autogenerated on Sun Oct 5 2014 22:29:52