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