#!/usr/bin/perl
#
# Do dynamic time warping (viterbi algorithm) on
# phoneme strings (list references).
#
# If not called within a require statement make a diff
# between two phoneme strings using a phoneme definition file
# like CGN_phonemes.pl:
# >perl dtw_phonemes.pl PhonemeDefinitionFile.pl <Reference> <Pattern>
#
# For example:
# >perl dtw_phonemes.pl CGN_phonemes.pl \
# 'kAbO+t@rs*e:t@n*ba:na:n@n' \
# 'k@bO+t@s*e:t@b@nAn@'
#
# Prints out a DIFF alignement linking Reference phonemen numbers
# with Pattern phoneme numbers (=) or insertion/deletions (|)
#
# Includes a printout of the split phoneme string on STDERR for confirmation
# Check this to see whether the phoneme defintion file is used correctly.
#
# Nice to find the word boundaries in connected speech when there is 
# a reference transcription. Use '*' to indicate a word boundary in the reference
# transcription.
# Please remove \" and \' quotes before processing.
#
# David Weenink, 2000
# Adapted by Rob van Son, 2001, 2004
#
#
#    Copyright (C) 2000, 2001  David Weenink, R.J.J.H. van Son
# 
#    This program 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.
# 
#    This program 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 this program; if not, write to the Free Software
#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#
#
# If you ALWAYS use CGN_phonemes.pl, uncomment the next line
#require 'CGN_phonemes.pl';

# Global vars for Viterbi.
# Costs associated with distance calculation for warping phoneme strings:
# vowel-consonant confusion e.g. i <-> l
$VCCost = 10;
# within Glides: l <-> r or Nasals: m <-> n or Vowels: iy <-> ih
$WithinCategoryCost = 3;
# between Categories, except vowel-consonant, e.g., Nasal-Glide: m <-> l
# !!! to do Fricative <-> Affricates more plausible than others
# !!! Plosive <-> Fricative more plausible than others
$BetweenRestCategoriesCost = 5;
#
# Substituting pauses is BAAAAAD
$PauseSubstitutionCost = 10;

# some directions and constants for the warp.
$DTW_DIAGONAL = 0; $DTW_HORIZONTAL = 1; $DTW_VERTICAL = -1;
$DTW_NODIRECTION = -2; $DTW_BIG = 1e38;

# Construct a distance matrix, i.e., the distance between ANY two phonemes
# in the reference and pattern strings.
sub distanceMatrix		# (\@ReferenceString, \@PatternString, \@ResultMatrix)
{
    my $reference = shift || return 0;
    my $pattern  = shift || return 0;
    my $distmatr  = shift || return 0;
    
    # Initialize
    my ($i, $j);
    
    # Initialize distance matrix
    @$distmatr = ();
    
    # Fill distance matric
    for ($i=0; $i <= $#$reference; $i++)
    {
	for ($j=0; $j <= $#$pattern; $j++)
	{
	    # The distance cost
	    $distmatr[$i][$j] = 
	    ($$pattern[$j] eq $$reference[$i]) ? 0 :	# No difference
	    sameCategory ($$pattern[$j], $$reference[$i]) ? $WithinCategoryCost : # Same category
	    (isOther ($$pattern[$j]) || isOther ($$reference[$i])) ? $PauseSubstitutionCost : # Pause			
	    (isVowel ($$pattern[$j]) || isVowel ($$reference[$i])) ? $VCCost : # Vowel to consonant
	    $BetweenRestCategoriesCost;	# Anything
	}
    };
    
    # Success
    return 1;
}   

# sub warpFind
#
# DTW of two phoneme strings (lists). 
# The optimal path and distance matrix are returned if references are entered.
#
# The return value is a DIFF like double list.
# The first list with the changes in the reference list, the second with the
# changes in the pattern list.
#
# input:
#	reference	:the reference phonemes as deduced from the dictionary (array)
#	pattern		:the realized phonemes (array)
# output:
#	path		:the optimal path  (Optional)
#	distmatr	:the distance matrix used for the warp (Optional)
#
# Result:
#	diff		: @$diff[0] = diff reference ([<>=|] <line>: <value>)
#			  @$diff[1] = diff pattern   ([<>=|] <line>: <value>)
#			  $diff[2]  = Total cost;

sub warpFind		# (\@reference, \@pattern [, \@path, \@distmatr]) -> @diff
{
    my $reference = shift || return 0;
    my $pattern   = shift || return 0;
    my $path      = shift || [];  # Optional
    my $distmatr  = shift || [];  # Optional
    
    # Initialize
    my ($i, $j); 
    my @delta = ();		# Summ costs
    my @psi = ();		# Store paths
    # Empty path
    @$path = ();
    
    # Calculate distance matrix
    &distanceMatrix ($reference, $pattern, $distmatr);
    
    # pattern and reference have same starting point (0,0) and ending point (last,last)
    $delta[0][0] = 0;
    # Fill delta and psi
    for ($i=1; $i <= $#$reference; $i++)
    {
	$delta[$i][0] = $DTW_BIG;
	$psi[$i][0] = $DTW_NODIRECTION;
    };
    
    #
    # Fill the cummulative cost, diff, and path matrices
    my $direction;
    my $minimum;
    for ($j=1; $j <= $#$pattern; $j++)
    {
	$delta[0][$j] = $delta[0][$j-1] + $distmatr[0][$j];
	$psi[0][$j] = $DTW_HORIZONTAL;
	for ($i=1; $i <= $#$reference; $i++)
	{
	    my $g;
	    $direction = $DTW_DIAGONAL;
	    $minimum = $delta[$i-1][$j-1] + 2 * $distmatr[$i][$j];
	    if (($g = $delta[$i][$j-1] + $distmatr[$i][$j]) < $minimum)
	    {
		$minimum = $g;
		$direction = $DTW_HORIZONTAL;
	    };
	    
	    if (($g = $delta[$i-1][$j] + $distmatr[$i][$j]) < $minimum)
	    {
		$minimum = $g;
		$direction = $DTW_VERTICAL;
	    };
	    
	    $psi[$i][$j] = $direction;
	    $delta[$i][$j] = $minimum;
	};	
    };
    
    # Get the best path and the diff of it
    my @diffref = ();
    my @diffpat = ();
    my $ipos = $#$reference;
    $minimum = $delta[$ipos][$#$pattern]; 
    # Start at the upper right corner of the path matrix
    my $index = $#$pattern + $#$reference + 1;
    $$path[$index][0] = $j = $#$pattern; #index 0= x; index 1= y
    $$path[$index][1] = $ipos;
    # Step down the matrix
    while ($j > 0)
    {
	if ($psi[$ipos][$j] == $DTW_DIAGONAL)
	{
	    if($$reference[$ipos] eq $$pattern[$j])
	    {
		push(@diffref, "= $ipos: $$reference[$ipos]");
		push(@diffpat, "= $j: $$pattern[$j]");
	    }
	    else
	    {
		push(@diffref, "< $ipos: $$reference[$ipos]");
		push(@diffpat, "> $j: $$pattern[$j]");
	    };
	    
	    $j--; $ipos--;
	}
	elsif ($psi[$ipos][$j] == $DTW_HORIZONTAL)
	{
	    push(@diffref, "| $ipos: $$reference[$ipos]");
	    push(@diffpat, "> $j: $$pattern[$j]");
	    
	    $j--;
	}
	elsif ($psi[$ipos][$j] == $DTW_VERTICAL)
	{
	    push(@diffref, "< $ipos: $$reference[$ipos]");
	    push(@diffpat, "| $j: $$pattern[$j]");
	    
	    $ipos--;
	}
	else
	{
	    die "warp: illegal path j = $j";
	};
	
	$$path[--$index][0] = $j;
	$$path[$index][1] = $ipos;
    };
    #
    # Ready, take the first item
    if($$reference[$ipos] eq $$pattern[$j])
    {
	push(@diffref, "= $ipos: $$reference[$ipos]");
	push(@diffpat, "= $j: $$pattern[$j]");
    }
    else
    {
	push(@diffref, "< $ipos: $$reference[$ipos]");
	push(@diffpat, "> $j: $$pattern[$j]");
    };
    
    # Reached the bottom, remove what is left on the lower axis
    while ($index-- > 0)
    {
	shift (@$path);
    };
    
    # Ready
    @diffref = reverse(@diffref);
    @diffpat = reverse(@diffpat);
    my @diff = ();
    $diff[0] = \@diffref;
    $diff[1] = \@diffpat;
    $diff[2] = $minimum;
    
    return @diff;
    
}

# Print the optimal path (for checking the results)
sub warpPrint
{
    my $reference = shift || return 0;
    my $pattern  = shift || return 0;
    my $path      = shift || return 0;
    
    # Initialize
    my ($i, $j, @line);
    local $, = '';
    local $\ = "\n";
    
    # Start printing
    print "\n"; 
    
    # 
    for ($i = $#$reference; $i >= 0; $i--)
    {
	# Fill lines with "-"
	for ($j=0; $j <= $#$pattern; $j++)
	{
	    $line[$j] = "-";
	};
	# Print a "*" for the real path
	for ($j=0; $j <= $#$path; $j++)
	{
	    $line[ $$path[$j][0] ] = "*" if $$path[$j][1] == $i;
	};
	# Print the reference phoneme
	print $$reference[$i], " " x (5-length($$reference[$i])),
	(" "), @line; 
    };
    
    # Print the pattern phonemes
    for ($i=0; $i <= 2; $i++)
    {
	for ($j=0; $j <= $#$pattern; $j++)
	{
	    $line[$j] = substr ($$pattern[$j], $i, 1) || " ";
	};
	
	print " " x 6, @line; 
    };
}

#
# Print the diff information 
sub diffPrint		# (\@diff)
{
    my $diff = shift || return 0;
    
    my $i;
    for($i=0; $i < scalar(@{$$diff[0]}); ++$i)
    {
	    print "$$diff[0][$i]\t$$diff[1][$i]\n";
    };
    
    return 1;
    
}


#
# If not called within a require, compare two strings and print the diff
# The first argument should be the CGN_phonemes.pl phoneme definition file
unless(caller())
{
    my $Dictionary = shift(@ARGV);
    require $Dictionary;
    my @Reference = splitPhonemes($ARGV[0]);
    my @Pattern = splitPhonemes($ARGV[1]);
    
# Print Phoneme strings for confirmation
print STDERR "@Reference\n@Pattern\n\n";

    my @Diff = ();
    my (@path, @distmatr) = ();
    
    @Diff = warpFind(\@Reference, \@Pattern, \@path, \@distmatr);

    diffPrint(\@Diff);
    
}

# Make require happy
1;
