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