#!/usr/bin/python # -*- coding: utf-8 -*- # -------------------------------------------------- # File Name: dtw.py # Purpose: # Creation Date: 10-03-2017 # Last Modified: Thu, Mar 23, 2017 8:31:16 PM # Author(s): Mike Stout # Copyright 2017 The Author(s) All Rights Reserved # Credits: # -------------------------------------------------- import sys from mlpy import dtw_subsequence, dtw_std #from cpu_sDTW import dtw_subsequence import numpy as np import mlpy import matplotlib.pyplot as plt import matplotlib.cm as cm def read(a): with open(a) as f: #s= f.readlines() s= f.read() return s m = int(sys.argv[1]) n = int(sys.argv[2]) f1 = sys.argv[3] f2 = sys.argv[4] gapChar = sys.argv[5] def zscore(xs): xs = map(float, map(ord, xs)) mean = np.mean(xs) var = np.var(xs) return (np.array(xs)-mean)/ var x = read(f1)[m:n] x_ = zscore(x) y = read(f2)[m:n] y_ = zscore(y) # dst,_,pos = dtw_subsequence(x,y) #print dst,pos #print np.array(x_) #print np.array(y_) dist, cost, path = dtw_std(x_ ,y_, dist_only=False) dist, cost, path = dtw_subsequence(x_ ,y_) #def dist(xs,ys): return np.array([[ abs(x-y) for x in xs ] for y in ys ]) #dst = dist(x_,y_) #print dst #exit() fig = plt.figure(1) ax = fig.add_subplot(111) plot1 = plt.imshow(cost.T, origin='lower', cmap=cm.gray, interpolation='nearest') plot2 = plt.plot(path[0], path[1], 'w') xlim = ax.set_xlim((-0.5, cost.shape[0]-0.5)) ylim = ax.set_ylim((-0.5, cost.shape[1]-0.5)) plt.show() def splits(n, f, xs): return [ f+":\t" + xs[i:i+n] for i in range(0, len(xs), n)] def fix(xs, i, path): if i>0 and path[i]==path[i-1]: return gapChar else: return xs[path[i]] def recover(xs, path): return ''.join([ fix(xs, i, path) for i in xrange(len(path))]).replace("\n"," ") n = 200 x__ = recover(x,path[0]) y__ = recover(y,path[1]) aa = splits(200, f1, x__) bb = splits(200, f2, y__) zz = [ "\n"+a+"\n"+b for a,b in zip(aa,bb) ] for z in zz: print z def mrg((x,y)): return x+y # str((x,y)) def merge(xs,ys): return ''.join(map(mrg, zip(xs,ys))) with open("tmp",'w') as f: f.write(merge(x__, y__))