/*********************************************************************
DTWkernel.cpp (version 0.1)
Computes the global alignment kernel
Copyright 2006 Marco Cuturi
The corresponding paper for this code is
M.C, J.-P. Vert, O. Birkenes, T. Matsui,
"A kernel for time series based on global alignments",
and can be found in a preprint form at http://arxiv.org/abs/cs.CV/0610033
This code is a modification of the LAKernel original code
by Jean-Philippe Vert and Hiroto Saigo (C) 2005. See the webpage of J.-P. Vert for
more details.
DTWkernel is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
Foobar is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with Foobar; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*********************************************************************/
#include
#include
#include
#include
#include
#include
#include
#include
#include "tnt.h" /* This bit is to use the TNT package for arrays and matrices */
/* Useful constants */
#define LOG0 -10000 /* log(0) */
#define LOGP(x,y) (((x)>(y))?(x)+log1p(exp((y)-(x))):(y)+log1p(exp((x)-(y))))
double DTWkernelcompute(TNT::Array2D seq1 ,TNT::Array2D seq2, double sigma)
/* Implementation of the global alignment kernel which generalizes the DTW algorithm */
/* seq1 is a first sequence represented as a matrix of real elements. Each line i corresponds to the vector of observations at time i. */
/* seq2 is the second sequence formatted in the same way.
/* NOTE: we use the TNT::Array2D implementation of matrices. The code can be easily modified if you use other matrix objects.
/* NOTE: seq1 and seq2 may have a different number of lines but need to have the same number of columns. */
/* sigma stands for the bandwidth of the Gaussian kernel */
{
register int
i,j,ii, /* loop indexes */
cur, old, /* to indicate the array to use (0 or 1) */
curpos, frompos1,frompos2,frompos3; /* position in an array */
double aux , aux2;
int cl; /* length of a column for the dynamic programming */
int nX,nY,dimvect;
nX=seq1.dim1();
nY=seq2.dim1();
dimvect=seq1.dim2(); /* note that we assume that seq1.dim2() and seq2.dim2() are the same.
double sum=0;
TNT::Array2D grammat (nX,nY,0.0);
double S_=1/(sigma*sigma);
/*********************************************************/
/* Computation of the Gram matrix with a Gaussian kernel */
/*********************************************************/
for (i=0;i logM(2*cl);
/************************************************/
/* First iteration : initialization of column 0 */
/************************************************/
/* The log=proabilities of each state are initialized for the first column (x=0,y=0..nY) */
for (j=0;j