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 }
00074 else {
00075 if ( x == m_parms[mean] ) value = 1.0;
00076 }
00077 return value * m_parms[norm];
00078 }
00079
00080
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 }
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