00001 /* 00002 * HippoPlot PowerLaw class implementation 00003 * 00004 * Copyright (C) 2000, 2003 The Board of Trustees of The Leland 00005 * Stanford Junior University. All Rights Reserved. 00006 * 00007 * $Id: PowerLaw.cxx,v 1.3 2003/07/23 01:10:36 pfkeb Exp $ 00008 * 00009 */ 00010 00011 #ifdef HAVE_CONFIG_H 00012 #include "config.h" 00013 #else 00014 #ifdef _MSC_VER 00015 #include "msdevstudio/MSconfig.h" 00016 #endif 00017 #endif 00018 00019 #include "PowerLaw.h" 00020 00021 #include "FunctionHelper.h" 00022 00023 #include <cmath> 00024 #include <cassert> 00025 00026 #ifdef ITERATOR_MEMBER_DEFECT 00027 using namespace std; 00028 #else 00029 using std::max; 00030 using std::vector; 00031 #endif 00032 00033 PowerLaw::PowerLaw ( ) 00034 { 00035 initialize (); 00036 } 00037 00038 PowerLaw::PowerLaw ( double prefactor, double index ) 00039 { 00040 initialize (); 00041 00042 m_parms[0] = prefactor; 00043 m_parms[1] = index; 00044 } 00045 00046 void PowerLaw::initialize () 00047 { 00048 m_name = "PowerLaw"; 00049 m_parm_names.push_back ( "Prefactor" ); 00050 m_parm_names.push_back ( "Index" ); 00051 00052 resize (); 00053 } 00054 00055 FunctionBase * PowerLaw::clone () const 00056 { 00057 return new PowerLaw ( *this ); 00058 } 00059 00060 double PowerLaw::operator () ( double x ) const 00061 { 00062 return m_parms[0]*pow(x, m_parms[1]); 00063 } 00064 00065 /* virtual */ 00066 void 00067 PowerLaw:: 00068 initialParameters ( const FunctionHelper * helper ) 00069 { 00070 double min_x = helper->minCoord (); 00071 double max_x = helper->maxCoord (); 00072 max_x = (min_x + max_x)/2.; 00073 00074 double min_y, max_y; 00075 try { 00076 min_y = helper->valueAt(min_x); 00077 max_y = helper->valueAt(max_x); 00078 if (min_y != 0 && max_y != 0) { 00079 m_parms[1] = log( max_y/min_y ) / log( max_x/min_x ); 00080 m_parms[0] = max_y/pow(max_x, m_parms[1]); 00081 return; 00082 } 00083 } catch (...) { 00084 // do nothing 00085 } 00086 00087 // default behavior 00088 min_y = max(helper->minValue(), 1.); 00089 max_y = helper->maxValue(); 00090 m_parms[1] = log( max_y/min_y ) / log( max_x/min_x ); 00091 m_parms[0] = max_y/pow(max_x, m_parms[1]); 00092 } 00093 00094 double PowerLaw::derivByParm ( int i, double x ) const 00095 { 00096 switch ( i ) { 00097 case 0 : 00098 return pow(x, m_parms[1]); 00099 break; 00100 00101 case 1 : 00102 return m_parms[0]*pow(x, m_parms[1])*log(x); 00103 break; 00104 00105 default: 00106 assert ( false ); 00107 break; 00108 } 00109 return 0.0; 00110 }