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 }
00080 else {
00081 if ( x == m_parms[mean] ) value = 1.0;
00082 }
00083 return value * m_parms[norm];
00084 }
00085
00086
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 }
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 }
00168