Marsyas  0.2
/home/gperciva/src/marsyas/src/marsyas/RBF.cpp
00001 /*
00002 ** Copyright (C) 1998-2006 George Tzanetakis <gtzan@cs.uvic.ca>
00003 **  
00004 ** This program is free software; you can redistribute it and/or modify
00005 ** it under the terms of the GNU General Public License as published by
00006 ** the Free Software Foundation; either version 2 of the License, or
00007 ** (at your option) any later version.
00008 ** 
00009 ** This program is distributed in the hope that it will be useful,
00010 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 ** GNU General Public License for more details.
00013 ** 
00014 ** You should have received a copy of the GNU General Public License
00015 ** along with this program; if not, write to the Free Software 
00016 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00017 */
00018 
00019 #include "common.h"
00020 #include "RBF.h"
00021 #include <cmath>
00022 
00023 using namespace std;
00024 using namespace Marsyas;
00025 
00026 //#define MTLB_DBG_LOG
00027 
00028 RBF::RBF(mrs_string name):MarSystem("RBF", name)
00029 {
00030     addControls();
00031 }
00032 
00033 RBF::RBF(const RBF& a) : MarSystem(a)
00034 {
00035     ctrl_RBFtype_ = getctrl("mrs_string/RBFtype");
00036     ctrl_Beta_ = getctrl("mrs_real/Beta");
00037     ctrl_symmetricIn_ = getctrl("mrs_bool/symmetricIn");
00038 }
00039 
00040 RBF::~RBF()
00041 {
00042 }
00043 
00044 MarSystem* 
00045 RBF::clone() const 
00046 {
00047     return new RBF(*this);
00048 }
00049 
00050 void 
00051 RBF::addControls()
00052 {
00053     addctrl("mrs_string/RBFtype", "Gaussian", ctrl_RBFtype_);
00054     ctrl_RBFtype_->setState(true);
00055 
00056     addctrl("mrs_real/Beta", 1.0, ctrl_Beta_);
00057 
00058     addctrl("mrs_bool/symmetricIn", false, ctrl_symmetricIn_);
00059 }
00060 
00061 mrs_real 
00062 RBF::GaussianRBF(const mrs_real val) const
00063 {
00064     //as defined in:
00065     //http://en.wikipedia.org/wiki/Radial_basis_function
00066     return exp(-1.0*val*val*ctrl_Beta_->to<mrs_real>());
00067 }
00068 
00069 mrs_real
00070 RBF::MultiquadraticRBF(const mrs_real val) const
00071 {
00072     //as defined in:
00073     //http://en.wikipedia.org/wiki/Radial_basis_function
00074     mrs_real beta = ctrl_Beta_->to<mrs_real>();
00075     return sqrt(val*val+beta*beta);
00076 }
00077 
00078 mrs_real
00079 RBF::ThinPlateSplineRBF(const mrs_real val) const
00080 {
00081     //as defined in:
00082     //http://en.wikipedia.org/wiki/Radial_basis_function
00083     return val*val*log(val);
00084 }
00085 
00086 void
00087 RBF::myUpdate(MarControlPtr sender)
00088 {
00089     (void) sender;
00090     RBFtype_ = ctrl_RBFtype_->to<mrs_string>();
00091     if(RBFtype_ == "Gaussian")
00092         RBFfunc_ = &RBF::GaussianRBF;
00093     else if(RBFtype_ == "Multiquadratic")
00094         RBFfunc_ = &RBF::MultiquadraticRBF;
00095     else if(RBFtype_ == "ThinPlateSpline")
00096         RBFfunc_ = &RBF::ThinPlateSplineRBF;
00097     else
00098     {
00099         RBFfunc_ = NULL;
00100         MRSWARN("RBF::myUpdate - unsupported RBF function: " + RBFtype_);
00101     }
00102 
00103     ctrl_onObservations_->setValue(ctrl_inObservations_, NOUPDATE);
00104     ctrl_onSamples_->setValue(ctrl_inSamples_, NOUPDATE);
00105     ctrl_osrate_->setValue(ctrl_israte_, NOUPDATE); //[?]
00106     ostringstream oss;
00107     mrs_string inObsNames = ctrl_inObsNames_->to<mrs_string>();
00108     for (mrs_natural i = 0; i < inObservations_; ++i)
00109     {
00110         mrs_string inObsName;
00111         mrs_string temp;
00112         inObsName = inObsNames.substr(0, inObsNames.find(","));
00113         temp = inObsNames.substr(inObsNames.find(",")+1, inObsNames.length());
00114         inObsNames = temp;
00115         oss << "RBF_" << RBFtype_ << "_" << inObsName << ",";
00116     }
00117     ctrl_onObsNames_->setValue(oss.str(), NOUPDATE);
00118 }
00119 
00120 void 
00121 RBF::myProcess(realvec& in, realvec& out)
00122 {
00123     mrs_natural o,t;
00124     mrs_real res;
00125 
00126     if(ctrl_symmetricIn_->isTrue())
00127     {
00128         mrs_natural endLoop = min(inSamples_, inObservations_); // just to be sure...
00129         MRSASSERT(in.getRows () >= endLoop);
00130         MRSASSERT(in.getCols () >= endLoop);
00131         for(t=0; t<endLoop; ++t)
00132         {
00133             for(o=0; o<=t; ++o)
00134             {
00135                 if(this->RBFfunc_)
00136                 {
00137                     res = (this->*RBFfunc_)(in(o,t));
00138                     //check for NaN and Inf results
00139                     if(res != res)
00140                     {
00141                         MRSERR("RBF::myProcess(): calculation of RBF(x) @" << prefix_ << " is returning NaN/Inf for x(" << o << "," << t << ") = " << in(o,t));
00142                     }
00143                     else if(res == 0)
00144                     {
00145                         MRSERR("RBF::myProcess(): calculation of RBF(x) @" << prefix_ << " is returning 0 for x(" << o <<"," << t << ") = " << in(o,t));
00146                     }
00147                     out(o,t) = res;
00148                 }
00149                 else
00150                     out(o,t) = in(o,t); //just bypass data, unmodified
00151 
00152                 //symmetry
00153                 out(t,o) = out(o,t);
00154             }
00155         }
00156     }
00157     else
00158     {
00159         for(t=0; t<inSamples_; ++t)
00160         {
00161             for(o=0; o<inObservations_; ++o)
00162             {
00163                 if(this->RBFfunc_)
00164                 {
00165                     res = (this->*RBFfunc_)(in(o,t));
00166                     //check for NaN and Inf results
00167                     if(res != res)
00168                     {
00169                         MRSERR("RBF::myProcess(): calculation of RBF(x) @" << prefix_ << " is returning NaN/Inf for x(" << o << "," << t << ") = " << in(o,t));
00170                     }
00171                     else if(res == 0)
00172                     {
00173                         MRSERR("RBF::myProcess(): calculation of RBF(x) @" << prefix_ << " is returning 0 for x(" << o <<"," << t << ") = " << in(o,t));
00174                     }
00175                     out(o,t) = res;
00176                 }
00177                 else
00178                     out(o,t) = in(o,t); //just bypass data, unmodified
00179             }
00180         }
00181     }
00182 #ifdef MARSYAS_MATLAB
00183 #ifdef MTLB_DBG_LOG
00184     MATLAB_PUT(name_, "name");
00185     MATLAB_PUT(out, "out");
00186     MATLAB_EVAL("if (length(out)>1) figure(1);imagesc(out);title(name);colorbar; end");
00187 #endif
00188 #endif
00189 }