// lsqbm.mata 1.0 27 AUG 2007
// Mata libary for sq.pkg
// Boyer-Moore Algorithm for -egen sqfirstpos()- and -sqallpos()-
// Magdalena Luniak, WZB
// Support: luniak@wzb.eu

version 9.1
    mata:
    mata clear
    mata set matastrict on


//compute the "jump"-table
real matrix computeJumpTable(real vector alphabet, real vector pattern){

	real matrix jumpTable	//result
	real scalar i			//loop variable over the alphabet
	real scalar j			//loop variable over the pattern
	real scalar element		//current element of the pattern
	real scalar jump			//length of the jump
	real scalar lengthAlphabet


	lengthAlphabet = length(alphabet)
	//initialisation of jumptable
	jumpTable = J(lengthAlphabet, 2, .)

	//loop over the patternvector
	for(i=1; i<=lengthAlphabet; i++){
		element=alphabet[i]
		jump=0
		//search for the most right occurance of the element in the pattern
		for(j=length(pattern); j>=1 ; j--){
			if (pattern[j]==element){
			break
			}
			else{
			jump++
			}
		} //end loop over the patternvector

	//Jump distance will be saved in the second row of the jumpTable
		jumpTable[i,1]=element
		if(jump == 0){	//if jump == 0 the algorithm does not terminate
		jumpTable[i,2]=1
		}
		else{
		jumpTable[i,2]=jump
		}
	} //end loop over the alphabet

	return(jumpTable)
}



//return the length of the jump for the given element of alphabet

real scalar getJump(real matrix jumpTable, real scalar element){

	real scalar i		//loop over jumptable
	real jump			//length of the jump

	for(i=1; i<=rows(jumpTable); i++){
		if(element == jumpTable[i,1]) {
			jump = jumpTable[i,2]
			break
		} //if
	} //end loop variable

	return(jump)
} // end function


//save the length of each sequence in a vector
//

real vector lengthOfEachSeq(real vector id){
	real scalar counter		//counts number of sequences
	real scalar i 			//loop over the id-vector
	real scalar maxLength		//Maximal length of vector
	real vector lengthsV		//Vector of lengths of each sequence
	real scalar index 		//index in lengths
	real scalar number		// number of sequences
	real vector result		// result


	maxLength = length(id)
	counter = 1
	number = 0
	lengthsV = J(maxLength,1, .)
	index = 1

	for(i=1; i<=maxLength-1; i++){
		if(id[i]==id[i+1]){
			counter ++
		}else{
			lengthsV[index,1]=counter
			number++
			counter = 1
			index++
		} //end else
	} //end for

	//the last and secondlast feld
	if(id[maxLength]==id[maxLength-1]){
		lengthsV[index,1]=counter
		number++
	}else{
		//the length of last sequence is 1
		lengthsV[index,1]=1
		number++
	}

	//compression of result
	result = J(number, 1, .)
	for(i=1; i<=number; i++){
		result[i]= lengthsV[i]
	}

	return(result)
} //end lengthOfEachSeq


//returns the vector of indexes of the first occurances of the pattern for all observations
//assertion: inputdata are correct
//Parameter:
//real vector data - sequences
//real string idvar - variable of seuences' id
//
real vector getFirst(real vector data, real vector idvar, real vector pattern, real vector alphabet){

	real vector first			//result
	real matrix jumpTable		//jump-table
	real scalar i			//loop over the observations
	real vector eachSeqLength	//vector of the length of each sequence
	real scalar thisSeqLength	//the length of the current sequence
	real scalar seq			//navigation through the sequence to the right
	real scalar seq2			//navigation through the sequence to the left
	real scalar	firstIndex		//index of first element of the current sequence in the data
	real scalar	lastIndex		//index of last element of the current sequence in the data
	real vector sequence		//current examined sequence from datamatrix
	real scalar pat			//navigation through the pattern
	real scalar position		//position of the pattern in "sequence"
	real scalar patLength		//length of the pattern
	real scalar element		//element of the sequence
	real scalar counter		//counter of matched elements
	real scalar jump			//distance to jump


	jumpTable = computeJumpTable(alphabet, pattern)
	eachSeqLength = lengthOfEachSeq(idvar)
	patLength = length(pattern)

	first = J(length(data), 1, .)


	//loop over the length of each sequence

	firstIndex = 0
	lastIndex= 0
	for(i=1; i<=length(eachSeqLength); i++){

		thisSeqLength= eachSeqLength[i]
		firstIndex = lastIndex +1
		lastIndex = lastIndex + thisSeqLength
		sequence = data[|firstIndex\lastIndex|]
		seq= 1
		position = 0
	//search untill not found and sequence not empty

	while(position == 0 & seq <= length(sequence)){

		element = sequence[seq]


		// the last element matches
		if(element == pattern[patLength] && seq>=patLength) {

			counter = 1
			if(patLength == 1){
				position=seq
				break
			}else{
				seq2=seq-1

				for (pat=patLength-1; pat>=1; pat--) {

					if (sequence[seq2] == pattern[pat]) {
						counter++

						seq2--

						// if all elements match the pattern
						if(counter == patLength){
							position=seq2+1

							break
						} //if
					}	//if
					else
					{
						jump=getJump(jumpTable, sequence[seq])
						seq = seq +jump
						break
					} //else

				} //for
			} //else
		} //if
		//the last element does not match
		else
		{
			jump=getJump(jumpTable, sequence[seq])
			seq = seq+jump
		} //else

	} //while

	real scalar j
	for(j=firstIndex; j<=lastIndex; j++){
	first[j]=position
	}

	} //end loop over the id

return(first)
} //end function

//returns the vector of of indexes of the first occurances of the pattern for all observations
//check the correctness of input data

void BMFirst(
   string scalar seqvar,
   string scalar patvar,
   string scalar alpvar,
	string scalar idvar) {

	real vector data                //Dataset
	real vector patter             //The real vector for the patterstring including missings
	real vector alphabe            //The real vector for the alphabet including missings
	real vector pattern             //The real vector for the patterstring
	real vector alphabet            //The real vector for the alphabet
	real vector id	            //The real vector for the ids
 	real vector result		//result
	real scalar patLength		//length of pattern
	real scalar i 			//loop over the pattern
	real scalar j 			//loop over the alphabet
	real scalar included		//boolean: element of the pattern occures in the alphabet
	real scalar element		//element of the pattern
	real vector eachSeqLength	//vector of the length of each sequence


	// View the Data
      st_view(data=.,.,tokens(seqvar))
	st_view(patter=.,.,tokens(patvar))
	st_view(alphabe=.,.,tokens(alpvar))
	st_view(id=.,.,tokens(idvar))

	//trim
	pattern = vectortrim(patter)
	alphabet = vectortrim(alphabe)

	patLength = length(pattern)

	//The pattern is too long
	eachSeqLength = lengthOfEachSeq(id)
	if(max(eachSeqLength) < patLength){
	    errprintf("the length of the pattern (%f) is longer than the longest sequence in the data(%f)", patLength,max(eachSeqLength) )
		exit()
	}

	//The pattern consists of some symbols that are not included in the alphabet

	//loop over the pattern
	for(i=1; i<=patLength; i++){
		element = pattern[i]
		included = 0

		//loop over the alphabet
		for(j=1; j<= length(alphabet); j++){

			if (element == alphabet[j]) {
			included = 1
			} //if
		} //end loop over the alphabet


		if (included == 0) {
			errprintf("The sign %f does not belong to the alphabet of sequences", element )
			exit()
		} // if

	} //end loop over the pattern


	// Store back to Stata
	result = getFirst(data, id, pattern, alphabet)
	st_addvar("long","_SQBMFirst")
        st_store(.,"_SQBMFirst",result)
} //end function



//removes missing values from the tail of the vectror
//and reduces the length of vector
	real vector vectortrim(real vector vect){
		real scalar i
		real scalar counter
		real vector result

		counter=0
	//Number of valid values at the beginning
		for(i=1; i<length(vect); i++){
			if(vect[i]!=.){
				counter++
			}else{
				break
			} //else
		} //for
		result= vect[|1\counter|]
	return(result)
	}

//**********************************************all*******************************************

//returns the vector of number of the all independent occurances of the pattern for all observations (pattern)
//assertion: inputdata are correct
//Parameter:
//real vector data - sequences
//real string idvar - variable of sequences' ids
//
real vector getAll( real vector data, real vector idvar, real vector pattern, real vector alphabet){

	real vector first			//result
	real matrix jumpTable		//jump-table
	real scalar i			//loop over the observations
	real vector eachSeqLength	//vector of the length of each sequence
	real scalar thisSeqLength	//the length of the current sequence
	real scalar seq			//navigation through the sequence to the right
	real scalar seq2			//navigation through the sequence to the left
	real scalar	firstIndex		//index of first element of the current sequence in the data
	real scalar	lastIndex		//index of last element of the current sequence in the data
	real vector sequence		//current examined sequence from datamatrix
	real scalar pat			//navigation through the pattern
	real scalar number		//number of the patternoccurance in "sequence"
	real scalar patLength		//length of the pattern
	real scalar element		//element of the sequence
	real scalar counter		//counter of matched elements
	real scalar jump			//distance to jump


	jumpTable = computeJumpTable(alphabet, pattern)
	eachSeqLength = lengthOfEachSeq(idvar)
	patLength = length(pattern)

	first = J(length(data), 1, .)


	//loop over the length of each sequence

	firstIndex = 0
	lastIndex= 0
	for(i=1; i<=length(eachSeqLength); i++){

		thisSeqLength= eachSeqLength[i]
		firstIndex = lastIndex +1
		lastIndex = lastIndex + thisSeqLength
		sequence = data[|firstIndex\lastIndex|]
		seq= patLength
		number = 0
	//search untill sequence is not empty

	while(seq <= length(sequence)){
			element = sequence[seq]

		// the last element matches
		if(element == pattern[patLength] && seq>=patLength) {

			counter = 1
			if(patLength == 1){
				seq++
				number++
			}else{
			seq2=seq-1

			for (pat=patLength-1; pat>=1; pat--) {
				//checks, whether other elements match the pattern
				if (sequence[seq2] == pattern[pat]) {
					counter++
					seq2--

					// if all elements match the pattern
					if(counter == patLength){
						number++
						seq = seq+patLength
					} //if
				}	//if
				else
				{
					jump=getJump(jumpTable, sequence[seq])
					seq = seq +jump
					break
				} //else

			} //for
		} //if
		}
		//the last element does not match
		else
		{
			jump=getJump(jumpTable, sequence[seq])
			seq = seq+jump
		} //else

	} //while

	real scalar j
	for(j=firstIndex; j<=lastIndex; j++){
	first[j]=number
	}

	} //end loop over the id

return(first)
} //end function


void BMAll(
   string scalar seqvar,
   string scalar patvar,
   string scalar alpvar,
	string scalar idvar) {

	real vector data                //Dataset
	real vector patter             //The real vector for the patterstring including missings
	real vector alphabe            //The real vector for the alphabet including missings
	real vector pattern             //The real vector for the patterstring
	real vector alphabet            //The real vector for the alphabet
	real vector id	            //The real vector for the ids
 	real vector result		//result
	real scalar patLength		//length of pattern
	real scalar i 			//loop over the pattern
	real scalar j 			//loop over the alphabet
	real scalar included		//boolean: element of the pattern occures in the alphabet
	real scalar element		//element of the pattern
	real vector eachSeqLength	//vector of the length of each sequence


	// View the Data
      st_view(data=.,.,tokens(seqvar))
	st_view(patter=.,.,tokens(patvar))
	st_view(alphabe=.,.,tokens(alpvar))
	st_view(id=.,.,tokens(idvar))

	//trim
	pattern = vectortrim(patter)
	alphabet = vectortrim(alphabe)

	patLength = length(pattern)

	//The pattern is too long
	eachSeqLength = lengthOfEachSeq(id)
	if(max(eachSeqLength) < patLength){
	    errprintf("the length of the pattern (%f) is longer than the longest sequence in the data(%f)", patLength,max(eachSeqLength) )
		exit()
	}

	//The pattern consists of some symbols that are not included in the alphabet

	//loop over the pattern
	for(i=1; i<=patLength; i++){
		element = pattern[i]
		included = 0

		//loop over the alphabet
		for(j=1; j<= length(alphabet); j++){

			if (element == alphabet[j]) {
			included = 1
			} //if
		} //end loop over the alphabet


		if (included == 0) {
			errprintf("The sign %f does not belong to the alphabet of sequences", element )
			exit()
		} // if

	} //end loop over the pattern


	// Store back to Stata
	result = getAll(data, id, pattern, alphabet)
	st_addvar("long","_SQBMAll")
      st_store(.,"_SQBMAll",result)
} //end function


// Compile into a libary
	mata mlib create lsqbm, replace
	mata mlib add ///
	  lsqbm computeJumpTable() getJump() getFirst() BMFirst() ///
	lengthOfEachSeq() vectortrim() getAll() BMAll()
	mata mlib index



end
exit