/*********************************************************************
    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 <iostream.h>
#include <stdlib.h>
#include <string>
#include <valarray>
#include <stdio.h>
#include <time.h>
#include <unistd.h>
#include <vector.h>
#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 <double> seq1 ,TNT::Array2D <double> 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 <double> grammat (nX,nY,0.0);
  double S_=1/(sigma*sigma);

  /*********************************************************/
  /* Computation of the Gram matrix with a Gaussian kernel */
  /*********************************************************/
  for (i=0;i<nX;i++) {
	for (j=0;j<nY;j++) {
		sum=0;
			for (ii=0;ii<dimvect;ii++) {
				sum+=(seq1[i][ii]-seq2[j][ii])*(seq1[i][ii]-seq2[j][ii]);
			}
			grammat[i][j]=exp(-sum*S_);
		}
	}
  }	
    
  /* Initialization of the arrays */
  /* Each array stores two successive columns of the (nX+1)x(nY+1) table used in dynamic programming */
  cl = nY+1;           /* each column stores the positions in the aaY sequence, plus a position at zero */
  vector < double > 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<cl;j++) {
    logM[j]=LOG0;
  }
  logM[0]=0;

  /* Update column order */
  cur = 1;      /* Indexes [0..cl-1] are used to process the next column */
  old = 0;      /* Indexes [cl..2*cl-1] were used for column 0 */


  double logT=LOG0; // LOG OF THE FINAL KERNEL VALUE.

  /************************************************/
  /* Next iterations : processing columns 1 .. nX */
  /************************************************/

  /* Main loop to vary the position in aaX : i=1..nX */
  for (i=1;i<=nX;i++) {
    	/* Special update for positions (i=1..nX,j=0) */
    	curpos = cur*cl;                  /* index of the state (i,0) */
    	logM[curpos] = LOG0; 
    	/* Secondary loop to vary the position in aaY : j=1..nY */
    	for (j=1;j<=nY;j++) {
			curpos = cur*cl + j;            /* index of the state (i,j) */
			frompos1 = old*cl + j;            /* index of the state (i-1,j) */
			frompos2 = cur*cl + j-1;          /* index of the state (i,j-1) */
			frompos3 = old*cl + j-1;          /* index of the state (i-1,j-1) */

			/* Doing the updates, in two steps
			aux= LOGP (logM[frompos1],logM[frompos2] );
			logM[curpos] = LOGP( aux , logM[frompos3] ) + grammat[i-1][j-1];
		}  /* end of j=1:nY loop */
    /* Update the culumn order */
		cur = 1-cur;
		old = 1-old;
	  }
  }
  
  /* Return the logarithm of the kernel */
  return logM[curpos];
}