/ Code author: Yunyue Zhu, 09/2002.
/ Joint work with Dennis Shasha.
/	DTW search using PAA

INF:1e63
/BaseDist:{[x;y]:a*a:x-y}
BaseDist:{[x;y]:_abs x-y}

//new version of DTW, without stack
//|x|=|y|,the warping width is k
LocalDTW2:{[x;y;k]
	Table::((#x);1+2*k)# 0
	Table[0;k+!k+1]::+\BaseDist[x[0];y[!k+1]]
	i:1
	n:#x
	do[n-1
		offset:0
		do[1+2*k
			j:i+offset-k
			if[(j>-1)&(j<n)
				Table[i;offset]::BaseDist[x[i];y[j]]+&/(getTable[i;offset-1;k];getTable[i-1;offset+1;k];getTable[i-1;offset;k])
			]
			offset+:1
		]
		i+:1
	]
	:Table[n-1;k]
}

//new version of DTW, without stack
//|x| and |y| are not nessesary the same,the warping width is delta
//this is the fastest one I can get from k
//possible improvement:
//	reduce space requirement from m*n to 2*m*n*delta
//	but this will make the program much less readable and the coversion overhead could exceed the saving space
LocalDTW:{[x;y;delta]
	n:#x;m:#y;k:_ delta*m;kk:_ delta*n
	fraction:m%n;width:1+2*k
	Table:(n;m)# INF
	Trace::(n;m)# ,(-1;-1)
	
	Table[0;!(k+1)] :+\BaseDist[x[0];    y[!k+1]]
	Table[!(kk+1);0]:+\BaseDist[x[!kk+1];y[0]   ]
	i:1;jfloat:0
	do[n-1
		jfloat+:fraction;jint:_ jfloat
		offset:-k
		do[width
			j:offset+jint
			if[(j>0)&(j<m)
				Table[i;j]:BaseDist[x[i];y[j]]+&/(Table[i;j-1];Table[i-1;j];Table[i-1;j-1])
				min:&/(Table[i;j-1];Table[i-1;j];Table[i-1;j-1])
				if[Table[i;j-1]  =min;Trace[i;j]::(i;j-1)]
				if[Table[i-1;j]  =min;Trace[i;j]::(i-1;j)]
				if[Table[i-1;j-1]=min;Trace[i;j]::(i-1;j-1)]
			]
			offset+:1
		]
		i+:1
	]
	:Table[n-1;m-1]
}

LowerEnvelop:{[x;k]
	l:#x
 	x:(k # INF),x,(k # INF)
 	:&/+x[(!l)+\:(!1+2*k)]
}

UpperEnvelop:{[x;k]
	l:#x
 	x:(k # -INF),x,(k # -INF)
 	:|/+x[(!l)+\:(!1+2*k)]
}

Dist_Env:{[ll;uu;y]
	//ll, uu and y are of equal length, the upper and lower envlope of x
	:(+/(((ll-y)*(uu-y))>0)*(BaseDist[ll;y]) & (BaseDist[uu;y]))
}

//translate from (`pitch;`duration) to time series vector
vectorize2:{[note]
	i:0;s:()
	do[#note[0]
		s,:note[1;i] # note[0;i]
		i+:1
	]
	:s
}

NORM:1% _sqrt 2

BuildIndex:{[DB;NumOfFeature]
	Index:()
	i:0
	do[#DB;
		db_seq:DB[i]
		Index,: ,WHW[db_seq;NumOfFeature]
		i+:1
	]
	:Index
}

QueryGlobalTW:{[DB;query]
	lq:#query
	minDistance:INF
	
	i:0;
	do[#DB
		dist:GlobalDTW[DB[i];query]
		dist%:|/(lq;#DB[i])
		if[(dist<minDistance)
			minDistance:dist;
			index:i
		]		
		i+:1
	]
	:index
}

QueryLocalTW:{[DB;query;delta]
	lq:#query
	minDistance:INF
	
	i:0;
	do[#DB
		dist:LocalDTW[query;DB[i];delta]
		dist%:|/(lq;#DB[i])
		if[(dist<minDistance)
			minDistance:dist;
			index:i
			//\s:("i:"),($i)," dist:",($dist)
		]		
		i+:1
	]
	:index
}
	
FilterQuery:{[DB;query;delta]
	count1:0
	ll :LowerEnvelop[query;_ delta*#query]
	uu :UpperEnvelop[query;_ delta*#query]
	lq:#query
	minDistance:INF
	
	i:0
	do[#DB
		db_seq:DB[i]
		dist:UTW_Env2[ll;uu;db_seq]
		if[dist<minDistance
			count1+:1
			dist:LocalDTW[query;db_seq;delta]
			dist%:|/(lq;#db_seq)	
			if[(dist<minDistance)
				minDistance:dist;
				index:i
				//\s:("i:"),($i)," dist:",($dist)
			]
		]
		i+:1
		\s:"i:",$i
		\s:"dist:",$dist
	]
	\s:("count1:"),($count1)
	:index
}

FilterTopQuery:{[DB_t;DB_r;query;delta;topk]
	count1:0;count2:0
	width :_ delta*N
	le:LowerEnvelop[query;width]
	ue:UpperEnvelop[query;width]
	reduce_le:(+/+(-1;k) # le)%k
	reduce_ue:(+/+(-1;k) # ue)%k
	
	toplist:topk # -1
	topdist:topk # INF
	
	i:0
	do[#DB_t
		dist:Dist_Env[reduce_le;reduce_ue;DB_r[i]]%nn		/Envelope with dimension reduction using my method :PAA
		if[dist<topdist[topk-1]
			count1+:1
			dist:Dist_Env[le;ue;DB_t[i]]%N
			if[dist<topdist[topk-1]
				count2+:1
				/dist:LocalDTW2[query;DB_t[i];width]%N
				dist:LocalDTW[query;DB_t[i];delta]%N
				if[(dist<topdist[topk-1])
					pos: topdist _bin dist
					toplist:toplist[!pos],i,toplist[pos+!topk-pos+1]
					topdist:topdist[!pos],dist, topdist[pos+!topk-pos+1]
					//\s:("i:"),($i)," dist:",($dist)
				]
			]
		]
		i+:1
	]
	\s:("count1:"),($count1),("count2:"),($count2)
	:(toplist;topdist)
}

FilterQueryWithFeature:{[DB;INDEX;N;query;delta]
	count1:count2:0
	query_feature:dwtWarpDelta[query;N;delta]
	ll :LowerEnvelop[query;_ delta*#query]
	uu :UpperEnvelop[query;_ delta*#query]
	lq:#query
	minDistance:INF
	
	i:0
	do[#INDEX
		index_seq:INDEX[i]
		dist:UTW_Env2[query_feature[0];query_feature[1];index_seq]
		if[dist<minDistance
			count1+:1
			db_seq:DB[i]
			dist:UTW_Env2[ll;uu;db_seq]
			if[dist<minDistance
				count2+:1
				dist:LocalDTW[query;db_seq;delta]
				dist%:|/(lq;#db_seq)	
				if[(dist<minDistance)
					minDistance:dist;
					index:i
				]
			]
		]
		i+:1
	]
	\s:("count1:"),($count1),(" count2:"),($count2)
	:index
}

Demo:{
	N:2
	db:1:"MusicNotes2.l"
	list:buildList[db]
	INDEX:BuildIndex[list;N]
	llist:#list
	hh:*1 _draw #list
	sublist:(!hh),(hh+1+!llist-hh+1)
	
	//\t h:  FilterQuery[list[sublist];list[hh];0.1]		
	//("ans:"),($h)
	//\t h: QueryLocalTW[list[sublist];list[hh];0.1]	
	//("ans:"),($h)
	//\t h:QueryGlobalTW[list[sublist];list[hh]]	 	
	//("ans:"),($h)
	//\t h:FilterQueryWithFeature[list[sublist];INDEX[sublist];N;list[hh];0.1]
	//("ans:"),($h)
}


test:{[hh;hhh]
	delta:0.1
	query:list[hh]
	ll :LowerEnvelop[query;_ delta*#query]
	uu :UpperEnvelop[query;_ delta*#query]
	lq:#query
	
	db_seq:list[hhh]
	\dist:UTW_Env2[ll;uu;db_seq]
	\dist:LocalDTW[query;db_seq;delta]
	\dist%:|/(lq;#db_seq)	
}