00001 00012 #ifdef _MSC_VER 00013 #include "msdevstudio/MSconfig.h" 00014 #endif 00015 00016 #include "Weibull.h" 00017 00018 #include "FunctionHelper.h" 00019 00020 #include <cmath> 00021 #include <cassert> 00022 #include <iostream> 00023 00024 #ifdef ITERATOR_MEMBER_DEFECT 00025 using namespace std; 00026 #else 00027 using std::vector; 00028 #endif 00029 00030 namespace hippodraw { 00031 00032 Weibull::Weibull ( ) 00033 { 00034 initialize (); 00035 } 00036 00037 Weibull::Weibull ( double prefactor, double scale, double shape ) 00038 { 00039 initialize (); 00040 00041 m_parms[0] = prefactor; 00042 m_parms[1] = scale; 00043 m_parms[2] = shape; 00044 } 00045 00046 void Weibull::initialize () 00047 { 00048 m_name = "Weibull"; 00049 m_parm_names.push_back ( "Prefactor" ); 00050 m_parm_names.push_back ( "Scale" ); 00051 m_parm_names.push_back ( "Shape" ); 00052 00053 resize (); 00054 } 00055 00056 FunctionBase * Weibull::clone () const 00057 { 00058 return new Weibull ( *this ); 00059 } 00060 00061 double Weibull::operator () ( double x ) const 00062 { 00063 return m_parms[0]*exp( -pow(x/m_parms[1],m_parms[2]) )*pow(x,m_parms[2] - 1.0); 00064 } 00065 00066 /* virtual */ 00067 void 00068 Weibull:: 00069 initialParameters ( const FunctionHelper * helper ) 00070 { 00071 double min_x = helper->minCoord (); 00072 double max_x = helper->maxCoord (); 00073 max_x = (min_x + max_x)/2.; 00074 00075 m_parms[2] = 1.0; 00076 00077 try { 00078 double min_y = helper->valueAt (min_x); 00079 double max_y = helper->valueAt (max_x); 00080 if (min_y != 0 && max_y != 0) { // success! 00081 m_parms[1] = ( min_x - max_x ) / log( max_y/min_y ); 00082 m_parms[0] = max_y / exp( -max_x/m_parms[1] ); 00083 return; 00084 } 00085 } catch (std::string &message) { 00086 std::cerr << message << std::endl; 00087 } 00088 00089 // All cleverness fails, so use default values.... 00090 m_parms[0] = 1.; 00091 m_parms[1] = 1.; 00092 } 00093 00094 double Weibull::derivByParm ( int i, double x ) const 00095 { 00096 switch ( i ) { 00097 case 0 : 00098 return operator()(x)/m_parms[0]; 00099 break; 00100 00101 case 1 : 00102 return operator()(x)*pow(x/m_parms[1],m_parms[2]) * m_parms[2]/m_parms[1]; 00103 break; 00104 00105 case 2 : 00106 return operator()(x)*(log(x) - pow(x/m_parms[1],m_parms[2]) * log(x/m_parms[1])); 00107 break; 00108 00109 default: 00110 assert ( false ); 00111 break; 00112 } 00113 return 0.0; 00114 } 00115 00116 } // namespace hippodraw 00117