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 using namespace hippodraw;
00033 
00034 Gaussian::Gaussian ( )
00035 {
00036   initialize ();
00037 }
00038 
00039 Gaussian::Gaussian ( double n, double m, double s )
00040 {
00041   initialize ();
00042   
00043   m_parms[norm] = n;
00044   m_parms[mean] = m;
00045   m_parms[sigma] = s;
00046 }
00047 
00048 void Gaussian::initialize ()
00049 {
00050   m_name = "Gaussian";
00051 
00052   m_parm_names.push_back ( "Norm" );
00053   m_parm_names.push_back ( "Mean" );
00054   m_parm_names.push_back ( "Sigma" );
00055 
00056   resize ();
00057 }
00058 
00059 FunctionBase * Gaussian::clone () const
00060 {
00061   return new Gaussian ( *this );
00062 }
00063 
00064 double Gaussian::operator () ( double x ) const
00065 {
00066   double value = 0.0;
00067   if ( m_parms[sigma] != 0.0 ) {
00068     double t = ( x - m_parms[mean] ) / m_parms[sigma];
00069     t = 0.5 * t*t;
00070     if ( fabs ( t ) < 50.0 ) {
00071       value = exp ( -t ) / ( 2.50662828 * m_parms[sigma] );
00072     }
00073   } // sigma == 0.
00074   else {
00075     if ( x == m_parms[mean] ) value = 1.0;
00076   }
00077   return value * m_parms[norm];
00078 }
00079 
00080 /* virtual */
00081 void 
00082 Gaussian::
00083 initialParameters ( const FunctionHelper * helper )
00084 {
00085   double min_x = helper->minCoord ();
00086   double max_x = helper->maxCoord ();
00087   int size = helper->size();
00088   double total = helper->getTotal ();
00089 
00090   m_parms[norm] = total * ( max_x - min_x ) / size;
00091   m_parms[mean] = helper->meanCoord ();
00092   m_parms[sigma] = helper->stdCoord ();
00093 }
00094 
00095 double Gaussian::derivByParm ( int i, double x ) const
00096 {
00097   switch ( i ) {
00098   case norm :
00099     return derivByNorm ( x );
00100     break;
00101 
00102   case mean :
00103     return derivByMean ( x );
00104     break;
00105 
00106   case sigma :
00107     return derivBySigma ( x );
00108     break;
00109 
00110   default :
00111     assert ( false );
00112     break;
00113   }
00114   return 0.0;
00115 }
00116 
00117 double Gaussian::derivByNorm ( double x ) const
00118 {
00119   if ( m_parms[sigma] != 0.0 ) {
00120     double t = ( x - m_parms[mean] ) / m_parms[sigma];
00121     t = 0.5 * t*t;
00122     if ( fabs ( t ) > 50.0 ) {
00123       return 0.0;
00124     }
00125     else {
00126       return exp ( -t ) / ( 2.50662828 * m_parms[sigma] );
00127     }
00128   } // sigma == 0.0
00129   else {
00130     if ( x != m_parms[mean] ) {
00131       return 0.0;
00132     } else {
00133       return 1.0;
00134     }
00135   }
00136 }
00137 
00138 double Gaussian::derivByMean ( double x ) const
00139 {
00140   double dx = x - m_parms[mean];
00141   if ( m_parms[sigma] != 0.0 ) {
00142     return m_parms[norm] * dx 
00143       * exp ( -dx*dx / ( 2.0 * m_parms[sigma] * m_parms[sigma] ) ) 
00144       / ( 2.50662828 * m_parms[sigma] * m_parms[sigma] * m_parms[sigma] );
00145   }
00146   else {
00147     if ( x != m_parms[mean] ) return 0.0;
00148     else return 1.0;
00149   }
00150 }
00151 
00152 double Gaussian::derivBySigma ( double x ) const
00153 {
00154   if ( m_parms[sigma] == 0.0 ) return 0.0;
00155   double dx = x - m_parms[mean];
00156   double p2 = m_parms[sigma] * m_parms[sigma];
00157   double ex = exp ( -dx*dx / ( 2.0 * p2 ) );
00158   return m_parms[norm] * ( dx*dx * ex / ( p2*p2) - ex / p2 ) / 2.50662828; 
00159 }
00160 

Generated for HippoDraw Class Library by doxygen