Gaussian.cxx

Go to the documentation of this file.
00001 
00012 #ifdef _MSC_VER
00013 #include "msdevstudio/MSconfig.h"
00014 #endif
00015 
00016 #include "Gaussian.h"
00017 
00018 #include "FunctionHelper.h"
00019 
00020 #include <cmath>
00021 #include <cassert>
00022 
00023 using std::distance;
00024 
00025 #ifdef ITERATOR_MEMBER_DEFECT
00026 using namespace std;
00027 #else
00028 using std::exp;
00029 using std::vector;
00030 #endif
00031 
00032 namespace {
00035   enum { norm = 0, mean = 1, sigma = 2 };
00036 }
00037 
00038 namespace hippodraw {
00039 
00040 Gaussian::Gaussian ( )
00041 {
00042   initialize ();
00043 }
00044 
00045 Gaussian::Gaussian ( double n, double m, double s )
00046 {
00047   initialize ();
00048   
00049   m_parms[norm] = n;
00050   m_parms[mean] = m;
00051   m_parms[sigma] = s;
00052 }
00053 
00054 void Gaussian::initialize ()
00055 {
00056   m_name = "Gaussian";
00057 
00058   m_parm_names.push_back ( "Norm" );
00059   m_parm_names.push_back ( "Mean" );
00060   m_parm_names.push_back ( "Sigma" );
00061 
00062   resize ();
00063 }
00064 
00065 FunctionBase * Gaussian::clone () const
00066 {
00067   return new Gaussian ( *this );
00068 }
00069 
00070 double Gaussian::operator () ( double x ) const
00071 {
00072   double value = 0.0;
00073   if ( m_parms[sigma] != 0.0 ) {
00074     double t = ( x - m_parms[mean] ) / m_parms[sigma];
00075     t = 0.5 * t*t;
00076     if ( fabs ( t ) < 50.0 ) {
00077       value = exp ( -t ) / ( 2.50662828 * m_parms[sigma] );
00078     }
00079   } // sigma == 0.
00080   else {
00081     if ( x == m_parms[mean] ) value = 1.0;
00082   }
00083   return value * m_parms[norm];
00084 }
00085 
00086 /* virtual */
00087 void 
00088 Gaussian::
00089 initialParameters ( const FunctionHelper * helper )
00090 {
00091   double min_x = helper->minCoord ();
00092   double max_x = helper->maxCoord ();
00093   int size = helper->size();
00094   double total = helper->getTotal ();
00095 
00096   m_parms[norm] = total * ( max_x - min_x ) / size;
00097   m_parms[mean] = helper->meanCoord ();
00098   m_parms[sigma] = helper->stdCoord ();
00099 }
00100 
00101 double Gaussian::derivByParm ( int i, double x ) const
00102 {
00103   switch ( i ) {
00104   case norm :
00105     return derivByNorm ( x );
00106     break;
00107 
00108   case mean :
00109     return derivByMean ( x );
00110     break;
00111 
00112   case sigma :
00113     return derivBySigma ( x );
00114     break;
00115 
00116   default :
00117     assert ( false );
00118     break;
00119   }
00120   return 0.0;
00121 }
00122 
00123 double Gaussian::derivByNorm ( double x ) const
00124 {
00125   if ( m_parms[sigma] != 0.0 ) {
00126     double t = ( x - m_parms[mean] ) / m_parms[sigma];
00127     t = 0.5 * t*t;
00128     if ( fabs ( t ) > 50.0 ) {
00129       return 0.0;
00130     }
00131     else {
00132       return exp ( -t ) / ( 2.50662828 * m_parms[sigma] );
00133     }
00134   } // sigma == 0.0
00135   else {
00136     if ( x != m_parms[mean] ) {
00137       return 0.0;
00138     } else {
00139       return 1.0;
00140     }
00141   }
00142 }
00143 
00144 double Gaussian::derivByMean ( double x ) const
00145 {
00146   double dx = x - m_parms[mean];
00147   if ( m_parms[sigma] != 0.0 ) {
00148     return m_parms[norm] * dx 
00149       * exp ( -dx*dx / ( 2.0 * m_parms[sigma] * m_parms[sigma] ) ) 
00150       / ( 2.50662828 * m_parms[sigma] * m_parms[sigma] * m_parms[sigma] );
00151   }
00152   else {
00153     if ( x != m_parms[mean] ) return 0.0;
00154     else return 1.0;
00155   }
00156 }
00157 
00158 double Gaussian::derivBySigma ( double x ) const
00159 {
00160   if ( m_parms[sigma] == 0.0 ) return 0.0;
00161   double dx = x - m_parms[mean];
00162   double p2 = m_parms[sigma] * m_parms[sigma];
00163   double ex = exp ( -dx*dx / ( 2.0 * p2 ) );
00164   return m_parms[norm] * ( dx*dx * ex / ( p2*p2) - ex / p2 ) / 2.50662828; 
00165 }
00166 
00167 } // namespace hippodraw
00168 

Generated for HippoDraw Class Library by doxygen