PLearn 0.1
DTWKernel.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // DTWKernel.cc
00004 //
00005 // Copyright (C) 2007 Xavier Saint-Mleux, ApSTAT Technologies inc.
00006 //
00007 // Redistribution and use in source and binary forms, with or without
00008 // modification, are permitted provided that the following conditions are met:
00009 //
00010 //  1. Redistributions of source code must retain the above copyright
00011 //     notice, this list of conditions and the following disclaimer.
00012 //
00013 //  2. Redistributions in binary form must reproduce the above copyright
00014 //     notice, this list of conditions and the following disclaimer in the
00015 //     documentation and/or other materials provided with the distribution.
00016 //
00017 //  3. The name of the authors may not be used to endorse or promote
00018 //     products derived from this software without specific prior written
00019 //     permission.
00020 //
00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031 //
00032 // This file is part of the PLearn library. For more information on the PLearn
00033 // library, go to the PLearn Web site at www.plearn.org
00034 
00035 // Authors: Xavier Saint-Mleux
00036 
00040 #include "DTWKernel.h"
00041 
00042 namespace PLearn {
00043 using namespace std;
00044 
00045 PLEARN_IMPLEMENT_OBJECT(
00046     DTWKernel,
00047     "Kernel for Dynamic Time Warping",
00048     "Kernel for Dynamic Time Warping\n"
00049     "(see sect.4.7 of Rabiner, L. and Juang, B. 1993 'Fundamentals of Speech Recognition'. Prentice-Hall, Inc.)\n"
00050     );
00051 
00053 // DTWKernel //
00055 DTWKernel::DTWKernel()
00056     :subvec_length(-1), 
00057      local_paths(), 
00058      distance_type("L2"), 
00059      max_time_deviation(-1)
00060 {
00061     /*
00062      * default local paths: 
00063      * (0,1) or (1,0) cost 1.
00064      * (1,1) costs 2.
00065      */
00066     local_paths.push_back(TVec<LocalStep>(1,make_pair(make_pair(0,1),1.)));
00067     local_paths.push_back(TVec<LocalStep>(1,make_pair(make_pair(1,0),1.)));
00068     local_paths.push_back(TVec<LocalStep>(1,make_pair(make_pair(1,1),2.)));
00069 }
00070 
00072 // declareOptions //
00074 void DTWKernel::declareOptions(OptionList& ol)
00075 {
00076     declareOption(ol, "subvec_length", &DTWKernel::subvec_length,
00077                   OptionBase::buildoption,
00078                   "length of each sub-vec within an example");
00079 
00080     declareOption(ol, "local_paths", &DTWKernel::local_paths,
00081                   OptionBase::buildoption,
00082                   "Specifies local path constraints and weights. "
00083                   "In each path, steps should be listed from "
00084                   "the ending point to the starting point."
00085                   "e.g. if a path can start at (0,0), go to (1,1) "
00086                   "and then end at (2,1), it should be listed as: "
00087                   "[((1,0),0.5), ((1,1),0.5)] where '0.5' are weights");
00088 
00089     declareOption(ol, "distance_type", &DTWKernel::distance_type,
00090                   OptionBase::buildoption,
00091                   "Name of the 'distance' function to use "
00092                   "when comparing features (sub-vecs).\n"
00093                   "one of: 'L2'  : sqrt{sum{(x-y)^2}}\n"
00094                   "        'L1'  : sum{|x-y|}\n"
00095                   "        'pow2': sum{(x-y)^2}\n");
00096 
00097     declareOption(ol, "max_time_deviation", &DTWKernel::max_time_deviation,
00098                   OptionBase::buildoption,
00099                   "Maximum allowed difference between i and j; "
00100                   "negative means no limit.");
00101 
00102     inherited::declareOptions(ol);
00103 }
00104 
00105 void DTWKernel::declareMethods(RemoteMethodMap& rmm)
00106 {
00107     rmm.inherited(inherited::_getRemoteMethodMap_());
00108 
00109     declareMethod(
00110         rmm, "dtw", &DTWKernel::remote_dtw,
00111         (BodyDoc("Calc. DTW on two feature vectors\n"),
00112          ArgDoc("x1","first vector"),
00113          ArgDoc("x2","second vector"),
00114          RetDoc ("dpath, dpoint, bptrs")));
00115 }
00116 
00117 
00118 
00120 // build //
00122 void DTWKernel::build()
00123 {
00124     // ### Nothing to add here, simply calls build_
00125     inherited::build();
00126     build_();
00127 }
00128 
00130 // build_ //
00132 void DTWKernel::build_()
00133 {
00134     if(subvec_length <= 0)
00135         PLERROR("In DTWKernel::build_() : "
00136                 "subvec_length must be specified before build "
00137                 "(%d is not a valid value).", subvec_length);
00138 
00139     //distance function
00140     if(distance_type == "L2")
00141         dist_fn= L2distance;
00142     else if(distance_type == "L1")
00143         dist_fn= L1distance;
00144     else if(distance_type == "pow2")
00145         dist_fn= powdistance;
00146     else
00147         PLERROR("In DTWKernel::build_() : "
00148                 "distance_type should be one of: " 
00149                 "'L2', 'L1' or 'pow2' "
00150                 "('%s' is not a valid value).", distance_type.c_str());
00151 
00152     //calc. valid region from local_paths (min/max slope)
00153     TVec<LocalPath>::iterator it;
00154     TVec<LocalPath>::iterator itbeg= local_paths.begin();
00155     TVec<LocalPath>::iterator itend= local_paths.end();
00156     TVec<LocalStep>::iterator jt;
00157     TVec<LocalStep>::iterator jtend;
00158     bool first= true;
00159     for(it= itbeg; it != itend; ++it)
00160     {
00161         int i= 0, j= 0;
00162         jtend= it->end();
00163         for(jt= it->begin(); jt != jtend; ++jt)
00164         {
00165             i+= jt->first.first;
00166             j+= jt->first.second;
00167         }
00168         real rij= i;
00169         rij/= j;
00170         if(rij < slope_ij_min || first)
00171             slope_ij_min= rij;
00172         real rji= j;
00173         rji/= i;
00174         if(rji < slope_ji_min || first)
00175             slope_ji_min= rji;
00176         first= false;
00177     }
00178 
00179 }
00180 
00182 //   dtw    //
00184 void DTWKernel::dtw(const Vec& x1, const Vec& x2) const
00185 {
00186     PLASSERT(x1.length() % subvec_length == 0);
00187     PLASSERT(x2.length() % subvec_length == 0);
00188 
00189     int n1= x1.length() / subvec_length;// n1 features
00190     TVec<Vec> subvecs1(n1);
00191     for(int i= 0; i < n1; ++i)
00192         subvecs1[i]= x1.subVec(i*subvec_length, subvec_length);
00193     int n2= x2.length() / subvec_length;// n2 features
00194     TVec<Vec> subvecs2(n2);
00195     for(int j= 0; j < n2; ++j)
00196         subvecs2[j]= x2.subVec(j*subvec_length, subvec_length);
00197 
00198     //init: calc. point-to-point distances
00199     // also pre-calc. bounds on j for each i
00200     int i,j,jmin,jmax;
00201     real jmin0, jmax0;
00202     dpoint.resize(n1,n2);
00203     jbounds.resize(n1,2);
00204     for(i= 0; i < n1; ++i)
00205     {
00206         if(slope_ij_min == 0.)
00207         {// avoid div by zero... especially if 0/0
00208             jmin0= 0.;
00209             jmax0= n2;
00210         }
00211         else
00212         {
00213             jmin0= static_cast<real>(i+1-n1)/slope_ij_min + n2 - 1;
00214             jmax0= static_cast<real>(i)/slope_ij_min;
00215         }
00216         jmin= static_cast<int>(ceil(max(jmin0, slope_ji_min * i)));
00217         jmax= static_cast<int>(
00218             floor(min(jmax0, slope_ji_min*(i-n1+1) + n2 - 1)));
00219         if(max_time_deviation >= 0)
00220         {// max abs. delta time (if >= 0)
00221             jmin= max(jmin, i - max_time_deviation);
00222             jmax= min(jmax, i + max_time_deviation);
00223         }
00224         jbounds(i,0)= jmin;
00225         jbounds(i,1)= jmax;
00227         /*
00228          * TODO: make sure this is OK for exotic local_paths
00229          *       (also see: option desc. for local_paths)
00230          */
00232         for(j= jmin; j <= jmax; ++j)
00233             dpoint(i,j)= dist_fn(subvecs1[i], subvecs2[j]);
00234     }
00235 
00236     //recurs: calc. path distances
00237     dpath.resize(n1,n2);
00238     bptrs.resize(n1,n2);
00239     dpath(0,0)= dpoint(0,0); //starting point
00240     
00241     real mn; //min. found at each step
00242     int ai, aj; //'actual' coords when following a local path
00243     pair<int,int> scoords; //coords for a step
00244     real dist; //'distance' to a given point
00245     TVec<LocalPath>::iterator it;
00246     TVec<LocalPath>::iterator itbeg= local_paths.begin();
00247     TVec<LocalPath>::iterator itend= local_paths.end();
00248     TVec<LocalStep>::iterator jt;
00249     TVec<LocalStep>::iterator jtend;
00250     bool path_ok; //is this path valid?
00251     bool some_path_ok; //is there any valid path to those coords?
00252 
00253     for(i= 0; i < n1; ++i)
00254     {
00255         jmin= jbounds(i,0);
00256         if(i==0)
00257             jmin= max(1,jmin);//skip (0,0), already calc.
00258         jmax= jbounds(i,1);
00259         for(j= jmin; j <= jmax; ++j)
00260         {
00261             some_path_ok= false;
00262             mn= REAL_MAX;
00263             for(it= itbeg; it != itend; ++it)
00264             {// for all local paths
00265                 path_ok= true;
00266                 ai= i; aj= j;
00267                 dist= 0.;
00268                 jtend= it->end();
00269                 for(jt= it->begin(); jt != jtend; ++jt)
00270                 {// for each step in the local path
00271                     dist+= dpoint(ai,aj) * jt->second;
00272                     scoords= jt->first;
00273                     ai-= scoords.first;
00274                     aj-= scoords.second;
00275                     if(ai < 0 || aj < 0)
00276                     {
00277                         path_ok= false;
00278                         break;
00279                     }
00280                 }
00281                 if(ai < 0 || aj < jbounds(ai,0) || aj > jbounds(ai,1))
00282                     path_ok= false;
00283                 if(path_ok)
00284                 {
00285                     dist+= dpath(ai,aj);//add dist. to beg. of path
00286                     if(dist < mn || !some_path_ok)
00287                     {
00288                         mn= dist;
00289                         bptrs(i,j)= make_pair(ai,aj);
00290                         some_path_ok= true;
00291                     }
00292                 }
00293             }
00294             dpath(i,j)= mn;//will be REAL_MAX if no path... but should not happen
00295             if(!some_path_ok && i==n1-1 && j==n2-1)
00296                 PLERROR("In DTWKernel::dtw : can't reach end of path! "
00297                         "Check your local_paths, they may be inconsistent.");
00298         }
00299     }
00300 }
00301 
00303 // remote dtw //
00305 tuple<Mat, Mat, TMat<pair<int, int> > > DTWKernel::remote_dtw(const Vec& x1, const Vec& x2) const
00306 {
00307     dtw(x1,x2);
00308     return make_tuple(dpath, dpoint, bptrs);
00309 
00310 }
00312 // evaluate //
00314 real DTWKernel::evaluate(const Vec& x1, const Vec& x2) const 
00315 {
00316     dtw(x1, x2);
00317     return dpath.lastElement();
00318 }
00319 
00321 // makeDeepCopyFromShallowCopy //
00323 void DTWKernel::makeDeepCopyFromShallowCopy(CopiesMap& copies)
00324 {
00325     inherited::makeDeepCopyFromShallowCopy(copies);
00326 
00327     deepCopyField(local_paths, copies);
00328     deepCopyField(dpoint, copies);
00329     deepCopyField(dpath, copies);
00330     deepCopyField(bptrs, copies);
00331     deepCopyField(jbounds, copies);
00332 }
00333 
00334 } // end of namespace PLearn
00335 
00336 
00337 /*
00338   Local Variables:
00339   mode:c++
00340   c-basic-offset:4
00341   c-file-style:"stroustrup"
00342   c-file-offsets:((innamespace . 0)(inline-open . 0))
00343   indent-tabs-mode:nil
00344   fill-column:79
00345   End:
00346 */
00347 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines