Jan Martinovič
email: jan.martinovic@vsb.cz

Pavla Dráždilová
email: pavla.drazdilova@vsb.cz

Kateřina Slaninová
email: katerina.slaninova@vsb.cz

Tomáš Kocyan
email: tomas.kocyan@vsb.cz

Václav Snášel
email: vaclav.snasel@vsb.cz

Left-Right Oscillate Algorithm for Community Detection used in E-learning system

Learning management systems are widely used as a support of distance learning. Recently, these systems successfully help in present education as well. Learning management systems store large amount of data based on the history of users' interactions with the system. Obtained information is commonly used for further course optimization, finding e-tutors in collaboration learning, analysis of students' behavior, or for other purposes. The partial goal of the paper is an analysis of students' behavior in a learning management system. Students' behavior is defined using selected methods from sequential and process mining with the focus to the reduction of large amount of extracted sequences. The main goal of the paper is description of our Left-Right Oscillate algorithm for community detection. The usage of this algorithm is presented on the extracted sequences from the learning management system. The core of this work is based on spectral ordering. Spectral ordering is the first part of an algorithm used to seek out communities within selected, evaluated networks. More precise designations for communities are then monitored using modularity.

1 Introduction E-learning is a method of education which utilizes a wide spectrum of technologies, mainly internet or computer-based, in the learning process. It is naturally related to distance learning, but nowadays is commonly used to support face-to-face learning as well. Learning management systems (LMS) provide effective maintenance of particular courses and facilitate communication within the student community and between educators and students [START_REF] Dráždilová | Computational intelligence methods for data analysis and mining of elearning activities[END_REF]. Such systems usually support the distribution of study materials to students, content building of courses, preparation of quizzes and assignments, discussions, or distance management of classes. In addition, these systems provide a number of collaborative learning tools such as forums, chats, news, file storage etc.

Regardless of LMS benefits, huge amount of recorded data in large collections makes often too difficult to manage them and to extract useful information from them.

To overcome this problem, some LMS offer basic reporting tools. However, in such large amount of information the outputs become quite obscure and unclear. In addition, they do not provide specific information of student activities while evaluating the structure and content of the courses and its effectiveness for the learning process [START_REF] Zorrilla | Web usage mining project for improving web-based learning sites[END_REF]. The most effective solution to this problem is to use data mining techniques [START_REF] Castro | Applying data mining techniques to elearning problems[END_REF].

The main goal of the paper is the description of our Left-Right Oscillate algorithm for community detection. The usage of this algorithm is presented on the extracted sequences from a learning management system. The core of this work is based on spectral ordering. Spectral ordering is the first part of an algorithm used to seek out communities within selected, evaluated networks. More precise designations for communities are then monitored using modularity.

The discovery and analysis of community structure in networks is a topic of considerable recent interest in sociology, physics, biology and other fields. Networks are very useful as a foundation for the mathematical representation of a variety of complex systems such as biological and social systems, the Internet, the world wide web, and many others [START_REF] Dorogovtsev | Evolution of Networks: From Biological Nets to the Internet and WWW[END_REF][START_REF] Newman | The structure and dynamics of networks[END_REF]. A common feature of many networks is community structure, the tendency for vertices to divide into groups, with dense connections within groups and only sparser connections between them [START_REF] Girvan | Community structure in social and biological networks[END_REF][START_REF] Newman | Detecting community structure in networks[END_REF].

Analysis of Students' Behavior

Several authors published contributions with relation to mining data from e-learning systems to extract knowledge that describe students' behavior. Among others we can mention for example [START_REF] Hershkovitz | Learning about online learning processes and students' motivation through web usage mining[END_REF], where authors investigated learning process of students by the analysis of web log files. A 'learnograms' were used to visualize students' behavior in this publication. Chen et al. [START_REF] Chen | E-learning behavior analysis based on fuzzy clustering[END_REF] used fuzzy clustering to analyze e-learning behavior of students. El-Hales [START_REF] El-Halees | Mining students data to analyze learning behavior: a case study[END_REF] used association rule mining, classification using decision trees, E-M clustering and outlier detection to describe students' behavior. Yang et al. [START_REF] Yang | Construction and application of the learning behavior analysis center based on open e-learning platform[END_REF] presented a framework for visualization of learning historical data, learning patterns and learning status of students using association rules mining. The agent technology and statistical analysis methods were applied on student e-learning behavior to evaluate findings within the context of behavior theory and behavioral science in [START_REF] Chen | The evaluation and analysis of student e-learning behaviour[END_REF].

Our subject of interest in this paper is student behavior in LMS, which is recorded in form of events and stored in the logs. Thus, we can define the student behavior with the terms of process mining which are used commonly in business sphere. Aalst et al. [START_REF] Van Der Aalst | Discovering social networks from event logs[END_REF][START_REF] Van Der Aalst | Process Mining: Discovery, Conformance and Enhancement of Business Processes[END_REF] defines event log as follows:

Let A be a set of activities (also referred as tasks) and U as set of performers (resources, persons). E = A ×U is the set of (possible) events (combinations of an activity and performer). For a given set A, A * is the set of all finite sequences over A. A finite sequence over A of length n is mapping σ =< a 1 , a 2 , . . . , a n >, where a i = σ (i) for 1 ≤ i ≤ n. C = E * is the set of possible event sequences. A simple event log is a multiset of traces over A.

The behavioral patterns are discovered using similarity of extracted sequences of activities performed in the system.

A sequence is an ordered list of elements, denoted < e 1 , e 2 , . . . , e l >. Given two sequences α =< a 1 , a 2 , . . . , a n > and β =< b 1 , b 2 , . . . , b m >. α is called a subsequence of β , denoted as α ⊆ β , if there exist integers 1 ≤ j 1 < j 2 < . . . < j n ≤ m such that a 1 = b j1 , a 2 = b j2 , . . . , a n = b jn . β is than a super sequence of α.

For finding the behavioral patterns, we need to use the methods for the sequence comparison. There are generally known several methods for the comparison of two or more categorical sequences. On the basis of our previous work [START_REF] Slaninová | Dynamic time warping in analysis of student behavioral patterns[END_REF] we have selected the algorithm, witch deals with the different lengths of sequences and with the possible error or distortion inside the sequence.

Time-warped longest common subsequence (T-WLCS) [START_REF] Guo | Time-Warped Longest Common Subsequence Algorithm for Music Retrieval[END_REF] is the method, which combines the advantages of two methods from pattern mining -The longest common subsequence (LCSS) [START_REF] Hirschberg | Algorithms for the longest common subsequence problem[END_REF] and Dynamic time warping (DTW) [START_REF] Müller | Information Retrieval for Music and Motion[END_REF]. LCSS allows us to find the longest common subsequence of two compared sequences. DTW is used for finding the optimal visualization of elements in two sequences to match them as much as possible. Then, selected T-WLCS method is able to compare sequences of various lengths, it takes into consideration the order of elements in the sequences, and it is immune to minor distortions inside of one of the compared sequences. Moreover, the method emphasizes recurrence of the elements in one of the compared sequences.

To obtain behavioral patterns of similar sequences in LMS, we have used our proposed new Left-Right Oscillate algorithm for community detection. The algorithm is based on spectral ordering (see Section 3).

Spectral Clustering and Ordering

Spectral clustering has become one of the most popular modern clustering algorithms in recent years. It is one of the graph theoretical clustering techniques and is simple to implement, can be solved efficiently by standard linear algebra methods, and very often outperforms traditional clustering algorithms such as the k-means or single linkage (hierarchical clustering). A comprehensive introduction to the mathematics involved in spectral graph theory is the textbook of Chung [START_REF] Chung | Spectral Graph Theory[END_REF]. Spectral clustering algorithm uses eigenvalues and eigenvectors of Laplacian of similarity matrix derived from the data set to find the clusters. A practical implementation of the clustering algorithm is presented in [START_REF] Cheng | On a recursive spectral algorithm for clustering from pairwise similarities[END_REF]. Recursive spectral clustering algorithm is used in [START_REF] Dasgupta | Spectral clustering by recursive partitioning[END_REF]. There Dasgupta et al. analyzed the second eigenvector technique of spectral partitioning on the planted partition random graph model, by constructing a recursive algorithm. A spectral clustering approach to finding communities in graphs was applied in [START_REF] White | A spectral clustering approach to finding communities in graphs[END_REF].

Ding and He showed in [START_REF] Ding | Linearized cluster assignment via spectral ordering[END_REF] that a linear ordering based on a distance sensitive objective has a continuous solution which is the eigenvector of the Laplacian. Their solution demonstrate close relationship between clustering and ordering. They proposed direct K-way cluster assignment method which transforms the problem to linearization the clustering assignment problem. The linearized assignment algorithm depends crucially on an algorithm for ordering objects based on pairwise similarity metric. The ordering is such that adjacent objects are similar while objects far away along the ordering are dissimilar. They showed that for such an ordering objective function the inverse index permutation has a continuous (relaxed) solution which is the eigenvector of the Laplacian of the similarity matrix.

Modularity -Quality of Detected Communities

To quantify the quality of the subdivisions we can use modularity [START_REF] Newman | Finding and evaluating community structure in networks[END_REF], defined as the fraction of links between the nodes in the same community minus their expected value in a corresponding random graph [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]. Networks with the high modularity have dense connections between the nodes within community, but sparse connections between the nodes in the different communities. Modularity is often used in optimization methods for detecting community structure in the networks [START_REF] Newman | Modularity and community structure in networks[END_REF]. The value of the modularity lies in the range -0.5, 1 . It is positive, if the number of edges within groups exceeds the number expected on the basis of chance.

For a weighted graph G we have a weight function w : E → R. It is for example function of the similarity between the nodes v i and v j . The weighted adjacency matrix of the graph is the matrix W = (w i j ) i, j = 1, . . . , n. Than the degree of a vertex v i ∈ V in weighted graph is defined as

d i = n ∑ j=1 w i j .
The weighted degree matrix D is defined as the diagonal matrix with the weighted degrees d 1 , . . . , d n on the diagonal.

In terms of the edge weights, modularity

Q(C 1 , . . . ,C k ) is defined over a specific clustering into k known clusters C 1 , . . . ,C k as Q(C 1 , . . . ,C k ) = k ∑ i=1 (e ii - k ∑ j=1,i = j e i j )
where e i j = ∑ (u,v)∈E,u∈C i ,v∈C j w(u, v) with each edge (u, v) ∈ E included at most once in the computation.

Left-Right Oscillate Algorithm for Community Detection

Upon completing our study of various modifications of algorithms for spectral clustering, we designed our own algorithm for detecting communities within complex networks. This algorithm utilizes spectral ordering where similar vertices are closer to indexes and less similar vertices are further from indexes. When determining the ordering, it is necessary to calculate the eigenvector of the second smallest eigenvalue of the matrix L = D -W . Since we have designed our algorithm for large amounts of data in a complex network, we used Lanczos method to calculate the Fiedler vector. Once the Fiedler vector was calculated, we detected appropriate gaps that divide the vertices of a graph into communities. As observed in the experiment, this type of separation into gaps leads to several badly-assigned subgraphs. This is due to the fact that the Fiedler vector is only linear ordered, as is revealed in our data collection. The Left-Right algorithm (see Algorithm 3.1) we have designed for incorporating small subgraphs into larger communities, gradually increases modularity in a given calculation.

Algorithm 1 Left-Right Algorithm for Community Detection

Input: similarity matrix W = w i, j for i = 1, . . . n and S c size of smallest communities. Output: communities C k , modularity of detected communities 1. Create Laplacian L = D -W using a matrix of similarity W of a connected graph G = (V, E,W ). 2. Calculate the Fiedler vector (the second eigenvector of Laplacian). 3. Reorder vertices according to Fiedler vector. 4. Calculate the sums of antidiagonals Asum i = ∑ ∀ j w i-j,i+ j a Asum (i+1/2) = ∑ ∀ j w i-j,i+ j+1 for all i = 1, . . . n and determine sum i = Asum(i -1/2)/4 + Asum(i)/2 + Asum(i + 1/2)/4. 5. Approximate the discrete function sum i by its spline and determine its first and second derivation. Then find all local minimums and maximums. 6. Assign maximum gaps that lie between two local maximums. Divide the set of vertices according to its gaps. Obtain subsets SS k ⊂ V , where k = 1, . . . K is the amount of subsets. 7. Detect a community using the Left-Right Oscillate assigning algorithm (see Algorithm 2). Spectral ordering minimizes the sum of weighted edges multiplied to the power of the difference in index nodes with the edge incidence. The calculation used for this equation is the given eigenvector of the second smallest eigenvalue (Fiedler vector) matrix L = D -W . A visualized Fiedler vector and ordered matrix similarity (in agreement with the Fiedler vector) reveals the creation of several natural clusters which is assigned by our algorithm.

For finding the Fiedler vector of Laplacian above a large, sparse and symmetric matrix representative of the evaluated network, we used Lanczos method to partially solve the eigenvalue problem. To determine the dimension of Krylov subspaces (for a more precisely calculated Lanczos method), we used modularity for determining the quality of a detected community. In [START_REF] Ding | Linearized cluster assignment via spectral ordering[END_REF], there is an example of a symmetric Laplacian L sym = D -1 2 LD -1 2 . Our experiment revealed, that this solution is only appropriate for some networks.

Algorithm 2 Left-Right-Oscillate Assigned Input: subsets SS k ⊂ V , where k = 1, . . . K, S c size of smallest communities. Output: communities C k , modularity of detected partitioning.

1. Find connected components C j for subset SS 1 , which are greater than the selected size |C j | ≥ S c . These components create communities. We add the rest of the vertices v i ∈ SS 1 -C j to the next subset of vertices SS 2 . 2. Find next connected components C j for every subset SS k k = 2, . . . K , which are greater than the selected size |C j | ≥ S c . These components create communities. Attempt to assign other vertices to the previous community, which was established in the previous step. If the vertex has no edge leading to the previous community than we add the vertex to the next subset of vertices SS k+1 . Continue repeating this method 2 until reach the end of ordered vertices. 3. Going through through all subsets of vertices, connected components are assigned to C j for j = 1, . . . J -1 and C J contain a set of connected components smaller than the selected size. 4. Employ the same approach going right-left without "oscillation". Begin with

C J-1 = C J-1 ∪ C J .
The next step for the algorithm is to order indexes of vertices v i ∈ V for all i = 1, . . . n in compliance with ordering using Fiedler vector values. Because we want to find communities that are easily detected in a visual representation when ordered by a similarity matrix, we must determine where one community in a linear order ends and the next begins (find two nodes that belong to various communities). For this reason, we have calculated the value of antidiagonal sums above an ordered the set of vertices that capture a cluster overlap in neighboring vertices v i . We define cluster crossing as the sum of a small fraction of the pairwise similarities. This is aided by linear ordering data points. The goal is to find the nodes that lay in the areas with fewer edges. These vertices lie close to locales with minimum function that are attached by approximation of a cluster overlap discrete function Sum i . We assigned this approximation using the spline function, allowing for easy calculations of both the first and second derivation, which are used to assign local extremes. Between the two local maximum extremes of this function, there lie two vertices. In this area, these vertices represent a maximum gap (the difference in their Fiedler vector value). This gap determines the border between two potential communities.

Using this method, we have found the natural amount of 'communities' above a given evaluated network. Since the precision with which the Fiedler vector is calculated is a determining factor, and since in some cases vertices are incorrectly assigned, the result is an irrelevant component. The benefit of using our algorithm lies within its ability to assign isolated nodes (or very small subgraphs with selected sizes) to the nearest, most suitable, connected component that creates the nucleus of a future community. Within our assignments, we gradually arrive at a set of vertices V separated by gaps in individual subsets V k . If the found set V k does not create a connected subgraph (G k = (V k , E), we determine all connected components in this subgraph. The maximum connected subgraph then creates the nucleus of this community and all subgraphs smaller than the selected size are moved to the right. We then attempt to reassign the subgraph to the next subset of vertices V k+1 . Due to the linear nature of spectral or-dering, it is presumable that subgraphs not yet assigned are reordered to the next subset of vertices. This means that we add the vertices of these subgraphs to the vertices of the next subset (that came into existence along gaps and creates a subgraph of the original graph with a set of vertices V k+1 ). Then we test the connectivity of subgraph G k+1 , which was expanded by the nodes from the previous, unassigned subgraph. We go through the entire, spectrally ordered set of graph vertices employing this method. At the end of this process, we have created the most relevant of components within which we assign small subgraphs that are not yet assigned. Then, we repeat this approach in the opposite direction -going from right to left -and we try to add vertices for inspection in a subgraph. We may then assign the vertices to a connected subgraph with adjacency to a vertex of a given subgraph.

Once we assign a subset of vertices using gaps, and once we have detected connected components from left to right and vice versa, we always calculate the modularity for the obtained separation of graphs into subgraphs. Our results have revealed that our Left-right method increases modularity. The resulting connected subgraphs then create the structure of communities in the graph, which is demonstrated on well known data collection Zachary karate club in Table 1.

Sequence Extraction in LMS Moodle

In this section is presented the extraction of students' behavioral patterns performed in the e-learning educational process. The analyzed data collections were stored in the Learning Management System (LMS) Moodle logs used to support e-learning education at Silesian University, Czech Republic. The logs consist of records of all events performed by Moodle users, such as communication in forums and chats, reading study materials or blogs, taking tests or quizzes etc. The users of this system are students, tutors, and administrators; the experiment was limited to the events performed only by students.

Let us define a set of students (users) U, set of courses C and term Activity a k ∈ A, where A = P × B is a combination of activity prefix p m ∈ P (e.g. course view, resource view, blog view, quiz attempt) and an action b n ∈ B, which describes detailed information of an activity prefix (concrete downloaded or viewed material, concrete test etc.). Event e j ∈ E then represents the activity performed by certain student u i ∈ U in LMS. On the basis of this definition, we have created a set S i of sequences s i j for the user u i , which represents the students' (users') paths (sessions) on the LMS website. Sequence s i j is defined as a sequence of activities, for example s i j =< a 1 j , a 2 j , . . . , a q j >, which is j-th sequence of the user u i .

The sequences were extracted likewise the user sessions on the web; the end of the sequences was identified by at least 30 minutes of inactivity, which is based on our previous experiments [START_REF] Dráždilová | Creation of students' activities from learning management system and their analysis[END_REF]. Similar conclusion was presented by Zorrilla et al. in [START_REF] Zorrilla | Web usage mining project for improving web-based learning sites[END_REF].

Using this method, we have obtained a set of all sequences S = ∪ ∀i S i , which consisted of large amount of different sequences s l performed in LMS Moodle. We have selected the course Microeconomy A as an example for the demonstration of proposed method. In Table 2 is presented detailed information about the selected course. The obtained set S of sequences consisted of large amount of different sequences, often very similar. Such large amount of information is hard to clearly visualize and to present in well arranged way. Moreover, the comparison of users based on their behavior is computationally expensive with such dimension. Therefore, we present in the article [START_REF] Slaninová | Dynamic time warping in analysis of student behavioral patterns[END_REF] the identification of significant behavioral patterns based on the sequence similarity, which allows us to reduce amount of extracted sequences.

We have used T-WLCS methods for the similarity measurement of sequences. The T-WLCS find the longest common subsequence α of compared sequences β x and β y , where α ⊆ β x ∧ α ⊆ β y , with relation to T-WLCS. Similarity was counted by the Equation 1.

Sim(β x , β y ) = (l(α) * h) 2 l(β x ) * l(β y ) , (1)
where l(α) is a length of the longest common subsequence α for sequences β x and β y ; l(β x ) and l(β y ) are analogically lengths of compared sequences β x and β y , and

h = Min(l(β x ), l(β y )) Max(l(β x ), l(β y )) (2)
On the basis of selected T-WLCS method for finding the similarity of sequences, we have constructed the similarity matrix for sequences (|S| × |S|) which can be represented using tools of graph theory. For the visualization of network was constructed weighted graph G(V, E), where weight w is defined as function w : E(G) → R, when w(e) > 0. Set V is represented by set of sequences S, weights w are evaluated by the similarity of sequences, see Equation 1, depending on selected method. In Table 3 is more detailed description of weighted graphs of sequences, where weight is defined by T-WLCS method for selected threshold θ (treshold for edges filtering -edges with smaller weights are removed). The number of nodes for each graph is 5908.

Reduction of Large Amount of Sequences by Left-Right Oscillate Algorithm

We described the procedure for extraction of sequencec from the LMS system in previous parts of the paper. We created the graphs of sequences by T-WLCS method. The examples in this section show how the graphs of sequences are divided to clusters by Left-Right Oscillate algorithm. We will use the concept "clusters" instead of "communities" in this part of the paper because we used Left-Right Oscillate algorithm for finding clusters of sequences. For the illustration, we have selected three different graphs from Table 3. These are the graphs created with parameter θ is greater than 0.2, 0.4 and 0.7. In each graph, there was identified the largest connected component, and on the basis on this component was created the new graph. This graph was partitioned by our new Left-Right Oscillate algorithm (maximum cycles of Lanczos algorithm inside Left-Right Oscillate algorithm was set to 1500). The sizes of these newly generated graphs are presented in Table 4.

Individual outputs and quality cuts of graphs after Left-Right Oscillate algorithm can be seen in Table 5, Table 5 andTable 7. In these tables we have column "Modularity Type" where row "Original" is without applied Left-Right Algorithm and "Left-Right" is row with information after Left-Right algorithm. Other columns in these tables are "Modularity" (see section 3.1), "Clusters" with amount of clusters after partitioning and columns with sizes of the top three communities (columns "Size of 1th", "Size of 2th", "Size of 3th").

It is apparent that the modularity is improved, if the parameter S c (size of smallest clusters -see Algorithm 3.1) of the Left-Right Oscillate algorithm is greater than 1.

Conclusion

In the paper we introduced the Left-Right Oscillate algorithm, which allows us to improve the results of community detection based on spectral ordering. We showed effect of parameter S c on the quality of clustering of sequences, which were extracted from the Moodle e-learning system. This allows us to better identify the same behavior of students in the online e-learning system. Modularity was used for measuring the quality of the distribution of sequences within clusters. In the future work we want to consider using the Left-Right Oscillate algorithm for hierarchical graph of sequences partitioning. Thanks to this we want to aim a more appropriate division of student's behavioral patterns.

Fig. 1 .

 1 Fig. 1. Similarity Matrix and Permuted Similarity Matrix (Natural Number of Communities is 7)

Table 1 .

 1 Modularity Before and After Left-Right Algorithm for Zachary Karate Club Before Left-right Commun. After Left Commun. After Left-right Commun.

	Laplacian	0.272	11	0.361	4	0.361	4
	Normalized-cut 0.342	7	0.311	4	0.311	4

Table 2 .

 2 Description of Log File for Course Microeconomy A

	Records	Students	Prefixes	Actions	Sequences
	65 012	807	67	951	8 854

Table 3 .

 3 Description of Sequence Graphs for T-WLCS Method

				T-WLCS
	θ Isolated Nodes Edges Avg. Degree Avg. Weighted Degree
	0.1	31	5577366 944.036	179.431
	0.2	143	1739042 294.354	88.883
	0.3	606	534648	90.496	38.735
	0.4	1200	271826	46.010	22.998
	0.5	2465	103028	17.439	10.430
	0.6	3781	29298	4.959	3.596
	0.7	5038	8080	1.368	1.269
	0.8	5517	5914	1.001	0.997
	0.9	5568	5788	0.980	0.980

Table 4 .

 4 Description of Selected Biggest Clusters of Sequence Graphs for T-WLCS Method

	T-WLCS
	θ Nodes Edges
	0.2 5763 1739040
	0.4 4639 271732
	0.7 142	1030

Table 7 .

 7 Partitioning of Sequence Graph for θ ≤ 0.7 S c = 1 Modularity Type Modularity Clusters Size of 1th Size of 2nd Size of 3rd Modularity Type Modularity Clusters Size of 1th Size of 2nd Size of 3rd

	Original	0.26457	63	15	15	14
	Left-Right	0.26457	63	15	15	14
			S c = 3
	Modularity Type Modularity Clusters Size of 1th Size of 2nd Size of 3rd
	Original	0.26457	63	15	15	14
	Left-Right	0.38731	12	23	22	21
			S c = 6
	Original	0.26457	63	15	15	14
	Left-Right	0.46304	7	45	15	14
			S c = 10
	Modularity Type Modularity Clusters Size of 1th Size of 2nd Size of 3rd
	Original	0.26457	63	15	15	14
	Left-Right	0.45416	4	58	42	27

Acknowledgment

This work was partially supported by SGS, VSB -Technical University of Ostrava, Czech Republic, under the grant No. SP2012/151 Large graph analysis and processing and by the European Regional Development Fund in the IT4Innovations Centre of Excellence project (CZ.1.05/1.1.00/02.0070).


Krzysztof Misztal
email: krzysztof.misztal@fis.agh.edu.pl

Jacek Tabor
email: tabor@ii.uj.edu.pl

Khalid Saeed
email: saeed@agh.edu.pl

A new algorithm for rotation detection in iris pattern recognition

Keywords: iris pattern recognition, rotation estimation, rotation recovery, Fourier descriptors

A new method for finding the rotation angle in iris images for biometric identification is presented in this paper. The proposed approach is based on Fourier descriptors analysis and algebraic properties of vector rotation in complex space.

Introduction

The iris is an important candidate for the source of the unique human population characteristics. It is very stable over time -the iris is determined during the first years of our lives and does not change until the death. Moreover, the study shows the iris structure is minimally dependent on our genes and allows to identify even identical twins. In addition, the modern iris biometrics systems are fast and have a high accuracy verification.

In light of the above, it seems that the iris is ideal for biometrics. However, the iris recognition process is complex and the construction of a complete system requires addressing a number of problems. The current approaches to iris pattern recognition include: the Gabor wavelet approach by Daugman [START_REF] Daugman | High confidence visual recognition of persons by a test of statistical independence[END_REF], the Laplacian parameter approach by Wildes et al. [START_REF] Wildes | Iris recognition: an emerging biometric technology[END_REF], zero-crossings of the wavelet transform at various resolution levels by Boles et al. [START_REF] Boles | A human identification technique using images of the iris and wavelet transform[END_REF], the Independent Component Analysis approach by Huang et al. [START_REF] Huang | An efficient iris recognition system Machine Learning and Cybernetics[END_REF], the texture analysis using multi-channel Gabor filtering and wavelet transform by Zhu et al. [START_REF] Zhu | Biometric personal identification based on iris patterns[END_REF], and many others ideas. Each method has its own advantages and disadvantages.

Let us focus our attention on the classical John Daugman algorithm [START_REF] Daugman | How iris recognition works[END_REF] -the most popular approach to individual identification by iris pattern recognition. and discuss how the algorithm copes with rotations. Let us start with the explanation and role of iris rotation angle. Fig. 1 present rotated and non-rotated iris images. When capturing iris images, most systems seek to achieve the same environmental conditions and the position of the face. Under some assumed restrictions the obtained images provide the highest reliability of identification methods. However, these methods are not comfortable for the users. The system proposed by Daugman detects the rotation angle by "brute force", i.e., by testing all accessible angles and selects the best of them. However, such a way, as we will show later is not optimal in practice. Our aim in this paper is to present alternative approach to this method. The method determines the rotation angle directly. Based on this knowledge we hope to develop a more efficient algorithm to identify individuals.

In the first part of our work we briefly describe the classical Daugman algorithm. In the next chapter we will focus on the problem of detecting the rotation angle of the iris image. Then we describe the rotation angle estimation algorithm using Fourier transform. Further considerations will include construction of our algorithm. The last part of the work contains numerical experiments that demonstrate how our rotation angle estimation algorithm is relatively insensitive to both the eye and the camera:

rotation, contrast modification, illumination level modification, blurring, sharpening. This makes it suitable to be used even in the first stage of preprocessing images for the algorithm of individual verification.

Iris code matching algorithms survey

In this section we describe the main stages in the iris pattern recognition algorithm proposed by John Daugman [START_REF] Daugman | How iris recognition works[END_REF].

After positive iris detection on an image and automatic segmentation we can match the iris region. To simplify, we assume we can describe the iris using only two circles -the arcs of the upper and lower eyelids are invisible on the surface of the iris. The next step in the iris pattern recognition algorithm performs the normalization process, which includes conversion of the iris annulus into the rectangular sized fixed image.

In general, images standardization is the initial processing procedure which allows their mutual comparison and further analysis. Comparing methods use predetermined patterns of the iris therefore the iris images could be compared with patterns from the database. In the develop of authors' method, we relied on the method of the normalization, called Daugman's Rubber Sheet Model. This Afterwards, we are able to introduce the feature encoding methods. In order to describe and recognize large individuals population only the significant features must be encoded. In Daugman's algorithm Gabor filters are used to obtain an optimal coding. As an outcome a phase sequence is obtained, 2048 bytes of data to describe the phase characteristics of the iris in a polar coordinate system -IrisCode. It is not affected by contrast, camera gain or illumination levels, so it is common in iris pattern recognition. Our method also has such properties.

For matching, the Hamming distance is chosen as a metric for recognition. The Hamming distance is calculated using only the bits generated from the true iris region, and this modified Hamming distance formula is given as

HD = (code A ⊕ code B) ∩ mask A ∩ mask B mask A ∩ mask B
where codes A and B are the two bit-wise templates to compare, mask A and mask B are the corresponding noise masks.

In theory, Hamming distance of two iris templates generated from the same iris equals 0.0. However, since the normalization is not perfect this will not occur in practice. We obtain some noise that goes undetected, and hence a difference will always appear.

Daugman's approach -brute force matching

The following method, suggested by Daugman, corrects the misalignment in the normalized iris pattern. In order to cope with rotational inconsistencies, when the Hamming distance of two templates is calculated, one template is shifted to one bit left and right. Then the Hamming distance values are calculated. This bit offset in the horizontal direction corresponds to the primary market area of the iris angle indicated by the Rubber Sheet Model. As for the iris distance we choose only the lowest, which corresponds to the best match between two templates. Fig. 3: The shifting process for one shift left and right.

The number of bits transferred during each shift by two times the number of filters used, since each filter will generate two bits of information from one pixel normalized region. The real number of changes necessary to normalize rotational inconsistencies will be determined by the maximum angle difference between two pictures of the same eye. One change is defined as one shift to the left, then one shift to the right. The example of shifting process is presented in Fig. 3.

In practice, when we configure the algorithm, we fix the maximum number of shifts to the left and to the right. This brute force rotation angle estimation by comparisons is made with each element in the database -this is the main disadvantage of this method.

Authors' Approach: Fourier Descriptor Extended Methodangle calculation

Discrete Fourier Transform (DTF). For convenience of the reader and to established notation we recall the basic definitions of DFT. DFT transforms (c 0 , . . . , c N -1 ) ∈ C n into a sequence (C 0 , . . . , C N -1 ) of complex numbers using the frequency domain representation:

C k = 1 √ N N -1 n=0 c n exp -i 2πnk N , 0 ≤ k ≤ N -1,
where c n = a n + ib n , for k = 0, . . . , N -1. This representation is widely used in data compression [START_REF] Salomon | Data compression: the complete reference[END_REF], partial differential equations, data mining, etc. It is easy to see that

C 0 = 1 √ N N -1 n=0 c n and C 1 = 1 √ N N -1 n=0 c n exp -i 2πn N . (1)
Remark 1. Fourier descriptors can be used to capture the main details of a boundary. This property is valuable because the coefficients usually keep the shape information. Thus they can be used as the basis for differentiating between distinct boundary shapes [START_REF] Gonzalez | Digital Image Processing[END_REF]. This allows recovering the rotation for images with ordinary boundaries (see Fig. 4). However, it is useless in the case of rotation of the circles. Moreover, it can be affected by image modification like contrast correction, gamma correction, blurring, etc.

The Fourier descriptor approach to rotation detection in iris pattern recognition. In this section we present how to recover the rotation angle using the modified idea of Fourier transform with vector angle calculation. The method from the above remark is not suitable in this situation, because the iris is represented as concentric circles and their rotations cannot be detected by boundary positions testing.

Assume that after using Daugman's Rubber Sheet Model we obtain two normalized iris I and I R for original and rotated images respectively. Images are assumed to be in gray-scale and each pixel in the image is coded by a number from the set {0, . . . , 255}. Next, we resize the images to obtain the same size, for example we set the new size to 8 × 256 (image height × width).

As an outcome we obtain two matrices

Ĩ = (c kn ) k,n ∈ M 8×256 (IR), ĨR = (c R kn ) k,n ∈ M 8×256 (IR).
Next, we construct the vectors of features describing Ĩ and ĨR by Fourier descriptors. For each row of these matrices we calculate the first Fourier descriptor according to equation ( 1) We put

C k = 1
x = (C 1 , . . . , C 8 ), x R = (C R 1 , . . . , C R 8 )
. Under the considerations of subsection about Fourier Descriptor Method, we get that those vectors are in the form x = e iϕ v, x R = w, for some v, w ∈ C 8 . Remark 2. The problem of finding the angle of rotation in the iris images can be formulated mathematically as follows.

Let non-trivial vectors v, w ∈ C 8 be fixed. We are looking for the minimum of the function f :

f v,w : [0, 2π] ϕ → e iϕ v -w 2 ∈ IR + .
These functions describe the rotation one of the vectors relatively to the second one. We want to reduce the vectors to the same one-dimensional linear subspace. This will reduce the difference between them.

After simple transformation we get

f v,w (ϕ) = v 2 + w 2 -(e iϕ v, w + e iϕ v, w ) = v 2 + w 2 -2 (e iϕ v, w ).
It is easy to see that the minimum is reached when the value of (e iϕ v, w ) is the largest. A trivial observation shows that

(e iϕ v, w ) ≤ |e iϕ v, w | = | v, w |, for ϕ ∈ [0, 2π].
ϕ = Arg v, w | v, w | . ( 2
)
Thus, with a description by a vector of features, we can find the angle by which the iris has been rotated.

Therefore, we can apply observation from previous section and by equation ( 2) we obtain

ϕ = Arg x, x R | x, x R | , ( 3
)
where ϕ is the desired angle. Summarizing, our approach to the problem of rotation angle estimation is stated as follows:

1. perform Daugman's Rubber Sheet Model standardization for I and I R (we use rubber sheets of 8 × 256); 2. calculate Fourier Transform for each row of rubber sheets; 3. select the first Fourier descriptor from each rubber's row and build the vectors x and x R ; 4. estimate angle of rotation using equation (3); 5. compare irises by Hamming distance using estimated angle.

The above steps link properties of the Fourier transform and the algebraic calculations. The computation is based on simple mathematical operations which make it easy for implementation.

Application in iris matching

In the last part of this work we present the comparison of classical Daugman's shifting method with our algorithm. Moreover, we make some experiments with noising and other "destruction" of information on iris images (ex. blur, gamma correction, sharpen) and examine the algorithm under such circumstances [START_REF] Tomeo-Reyes | Input images in iris recognition systems: A case study[END_REF].

We prepare implementation of Daugman's iris recognition algorithm based on [START_REF] Daugman | How iris recognition works[END_REF] and [START_REF] Daugman | High confidence visual recognition of persons by a test of statistical independence[END_REF]. Besides we extend this application by adding the implementation of our method.

Comparison between classical and authors' method in rotation

angle detection.

Firstly, we compare the efficiency of shifting method with our method. Table 1 contains comparison between the classical method and ours. We evaluate IrisCode for rotated (for angle 0, . . . , 9 degrees) and non-rotated images. Then we run shifting method to see what angle (column Rec.) was chosen and what the value of Hamming distance (column HD) is. In this case IrisCode evaluated for rotated iris image was compared with shifted original IrisCode in maximum 10 shifts. Next, we use our method which estimates rotation angle (column Rec.), and then the estimated angle was rounded to the nearest integer and used to generate shifted IrisCode. The obtained distance is presented in column HD.

The experiment results show that our method gives the same angle as the brute force method, with the advantage that calculations were performed faster.

Sensitivity to image modifications

During the second experiment we want to check if our method is insensitive to image destruction. Fig. 5 presents images used in this experiment. But before we present experiment details let us put some notations.

Let I denote image of size w × h pixels (width × height) -it is a matrix, so I ∈ {0, . . . , 255} w×h for gray-scale image. Then I(x, y) for x ∈ {1, . . . , w}, y ∈ {1, . . . , h} denotes colors value at position (x, y) on image. By I N we denote the new image (output image) created by modifying the input image I.

We can now proceed in explaining the original image (Fig. 5a) modification descriptions. We use the following methods: blur (Fig. 5e) -it is a convolution filter with mask

B = 1 16   1 2 1 2 4 2 1 2 1   ,
used to reduce noise and details information on image. In this case linear filtering with a mask of size 3 × 3 is given by the expression

I N (x, y) = 1 s=-1 1 t=-1
B(s, t)I(x + s, y + t).

To generate full filtered image this equation should be applied for x ∈ {1, . . . , w} and y ∈ {1, . . . , h} (see [START_REF] Gonzalez | Digital Image Processing[END_REF]). sharpen (Fig. 5f) -it is a convolution filter with a mask

  -1 -1 -1 -1 9 -1 -1 -1 -1   .
Sharpens the image making it look crisper and making the edges in the image more distinct. Table 2 summarizes the values of estimated angles with error level obtained by our method. Table 3 contains angles recovered by classical Daugman method. The error level in both examples suggests that our method is not affected by the proposed image factors.

Conclusions

The iris identification and verification process is complex and the construction of a complete system requires addressing a number of problems. One of the most important and most difficult ones is that of detecting the rotation angle of the iris image. The authors have introduced a new method to solve this problem based on simple mathematical operations which have made it easy for implementation. Fourier transform and some algebraic calculations are utilized to work out an algorithm of high success rate whose performance is presented in this work. It has described and shown the rotation angle estimation using Fourier transform. The numerical experiments have demonstrated how the authors algorithm is relatively insensitive to both the eye and the camera relative to rotation, contrast changes, illumination level and other factors that usually affect the image under processing.

Fig. 1 :

 1 Fig. 1: Iris images: original (non-rotated) and rotated (23 degree).

Fig. 2 :

 2 Fig. 2: Daugman's Rubber Sheet Model

Fig. 4 :

 4 Fig. 4: Example of two shapes with digital boundary. The rotation of the second one does not allow to estimate angle of rotation.

,

  for k = 1, . . . , 8.

  Consequently, we get e iϕ = v, w /| v, w |, since z • z = |z| 2 for z ∈ C and |e iϕ | = 1. Thus the minimum of the function f is reached for

Fig. 5 :

 5 Fig. 5: Iris image modified by different methods: 5a -original image; 5b -brightness correction (+50); 5c -contrast correction (+30); 5d -gamma correction (+2.5); 5e -blur; 5f -sharpen.

Table 1 :

 1 Comparison of angle detection in classical methods with shifts and authors' method

		Classical	Authors'
	Angle Rec.	HD	Rec.	HD
	0	0	0	0.0	0
	1	0	0	0.0	0
	2	1	0	1.4	0
	3	3	0	2.8	0
	4	3	0	2.8	0
	5	4	0	4.2	0
	6	6	0	5.6	0
	7	6	0	5.6	0
	8	7	0	7.0	0
	9	8	0	8.4	0

Table 2 :

 2 Angle detection in different image modifications (angle is given in degrees).

				Angle
	Method	5	14	23	47	76
		Rec. Error Rec. Error Rec. Error Rec. Error Rec. Error
	None	4.22 0.78 12.66 1.34 22.50 0.50 46.41 0.59 75.94 0.06
	Brightness 4.22 0.78 12.66 1.34 22.51 0.49 46.41 0.59 75.94 0.06
	Contrast 3.99 1.01 12.43 1.57 22.27 0.73 46.18 0.82 75.71 0.29
	Gamma 3.98 1.02 12.41 1.59 22.26 0.74 46.16 0.84 75.70 0.30
	Blur	3.69 1.31 11.99 2.01 21.87 1.73 45.69 1.31 75.13 0.87
	Sharpen 3.66 1.34 13.72 0.28 23.10 0.10 47.92 0.92 77.83 1.83

Table 3 :

 3 Angle detection in different image modifications (angle is given in degrees) obtained by classical shifting.

						Angle
	Method	5		14			23	47		76
		Rec. Error Rec. Error Rec. Error Rec. Error Rec. Error
	None	5	0	14	0	23	0	47	0	76	0
	Brightness 0	5	10	4	0	23	38	9	8	68
	Contrast	0	5	10	4	0	23	38	9	8	68
	Gamma	0	5	11	3	1	22	38	9	8	68
	Blur	0	5	13	2	3	20	60 13	30 46
	Sharpen	0	5	10	4	0	23	38	9	8

Acknowledgment

The work is supported by AGH University of Science and Technology, Krakow (grant no. 11.11-220-01).


Mariusz Rybnik
email: mariuszrybnik@wp.pl

Wladyslaw Homenda
email: homenda@mini.pw.edu.pl

Tuning of a Knowledge-driven Harmonization Model for Tonal Music

Keywords: harmonization, data analysis, expert system, musical work, supervised learning, tonal music

The paper presents and discusses direct and indirect tuning of a knowledge-driven harmonization model for tonal music. Automatic harmonization is a data analysis problem: an algorithm processes a music notation document and generates specific meta-data (harmonic functions). The proposed model could be seen as an Expert System with manually selected weights, based largely on the music theory. It emphasizes universality -a possibility of obtaining varied but controllable harmonies. It is directly tunable by changing the internal parameters of harmonization mechanisms, as well as an importance weight corresponding to each mechanism. The authors propose also indirect model tuning, using supervised learning with a preselected set of examples. Indirect tuning algorithms are evaluated experimentally and discussed. The proposed harmonization model is prone both to direct (expert-based) and indirect (data-driven) modifications, what allows for a mixed learning and relatively easy interpretation of internal knowledge.

Introduction

Harmony is an important element of tonal music, it defines a vertical relation between notes [START_REF] Sikorski | Harmony[END_REF], and by definition is opposed to melody -a horizontal succession of notes in a specific voice. In fact, however, the harmonic relations of the leading melody (regarding mono-or homophony) are largely depending on the horizontal succession of notes and melodic intervals between them. Harmonic passages that follow harmonic chords (with smaller or larger deviations) can be frequently detected in a melody. Similar, but obviously much stronger harmonic relations are to be found in an accompaniment, where the melody is less (or even not at all) important, as the main goal is to define a background for the leading melody. The obvious exception to these assumptions are polyphonic musical works. They tend to cultivate two or more independent voices, that compete for attention, but also have to cooperate harmonically, usually indicating strong harmonic relations.

Automatic harmonization can be seen as a problem crossing two different areas: 1) theoretical music knowledge and 2) formal mathematical computations -presented in form of algorithms or slightly less formal Artificial Intelligence solutions. From the technical point of view it could be seen as analyzing of a music notation document and producing meta-data (harmonic functions). Over the years, many various approaches and techniques were used to solve the problem (or similar). The most popular paradigms are: Expert Systems [START_REF] Cope | An expert system for computer-assisted music composition[END_REF][3], Neural Networks [START_REF] Hild | Harmonet: A neural net for harmonizing chorals in the style of J.S. Bach[END_REF], Constraints and probabilistic approaches [START_REF] Pachet | Musical harmonization with constraints: A survey[END_REF] [6] [START_REF] Paiement | Probabilistic Melodic Harmonization[END_REF], evolutionary algorithms [START_REF] Prisco | An evolutionary music composer algorithm for bass harmonization[END_REF]. Due to excessive complication and uncertainty, related to such sensitive subject as music (feelings and sensitivity still tend to evade scientific approaches), no approach seems to fully explore and describe the subject. In fact, some of the approaches focus on a particular style of musical works, as for example typically Baroque pieces of J. S. Bach in [START_REF] Hild | Harmonet: A neural net for harmonizing chorals in the style of J.S. Bach[END_REF]. This allows for relatively narrow specialization and therefore eases the algorithmic description of the problem.

Frequent limitations for above-mentioned approaches seem to be: lack of universality (results limited to a specific area), an unpredictability, a need for a large learning examples database and an extensive amount of calculations. The authors' approach is mostly based on an advanced theoretical music knowledge, especially in the area of harmonization, and is aimed at solving these disadvantages. The authors propose a model that may be customized at various levels and does not need any learning examples, as they are replaced by an expert knowledge of harmony, incorporated in model. It may be defined initially as Expert System with manually selected weights, however the authors also propose a scheme for tuning the general harmonization model. With enough knowledge of harmony one can also directly modify the model using internal parameters.

In this paper the authors concentrate on tuning procedures of their harmonization model introduced in [START_REF] Rybnik | Knowledge-driven Harmonization Model for Tonal Music[END_REF] and extended in [START_REF] Rybnik | Extension of Knowledge-driven Harmonization Model for Tonal Music[END_REF]. The paper is organized into five sections; Section 1 being an introduction. Section 2 presents basic concepts concerning harmony and harmonization in tonal music. Section 3 describes the harmonization model and explains mechanisms used. Section 4 discuses various tuning possibilities of a proposed model and presents a sample of experimental results. Finally, Section 5 concludes the paper, discuses configuration, universality properties and suggests future work.

On harmony in tonal music

The process of creating accompaniment to a homophony (lone melody -a singlelayered progression of sounds in time) can, in practice, be divided into two phases. The first is harmonization (determination of harmonic functions and chords corresponding to them), the second phase is creation of accompaniment using previously determined chords. It is important to stress, that in vast number of cases there is no single, ultimate harmonization for a given melody. There are many possibilities, from the most simple and obvious to the highly complicated and non-trivial ones (e.g. improvised music, jazz). The decision depends mostly on a style of music. It also relies on capabilities of the available instrument/orchestration, and the abilities of performers/harmonizers. There are also cases when no additional information is needed to perform the music, except for the leading melody and harmonic functions. A common example is an improvised à vista accompaniment to songs played on a guitar or a piano.

The tonal harmony

For the purposes of this paper the authors have focused on a tonal system with two basic scales: major and minor (natural minor as well as harmonic and melodic modifications). The authors' considerations and experiments are based on seven diatonic harmonic functions (built on 1 st , 2 nd , 3 rd , 4 th , 5 th , 6 th and 7 th grade of major/minor scales) with common modifications: adding the seventh, ninth, and sixth (all three minor and major).

The harmonic relations in tonal system tend to be well defined in terms of consonants (chords 'pleasant' to human ear, stable) and dissonances (chord slightly 'unpleasant', unstable, introducing a tension that needs to be resolved to consonants). The authors have omitted alterations and higher intervals (alike eleventh, thirteenth), that may be considered as dissonances meant to resolve, rather than intrinsic harmony.

Proposed harmonization model

The proposed harmonization model is based on several mechanisms that closely follow music theory:

1. Particular note can have various harmonic importance (based on the note's relative length, the notes placement in measure, surrounding notes, volume, etc.); 2. Each note excites (fits to) several harmonic functions, based on pitch and function components taken into consideration (components higher than 9 th are usually very rare and considering them is very difficult); 3. Some harmonic functions are more likely to occur than others (the simplest example being Tonic -a base and consolation for vast majority of tonal music, therefore usually occurring most frequently). It is preferable to prioritize the commonly used functions; 4. Some specific successions of harmonic functions are more or less probable, therefore it is possible to prioritize the more likeable (frequent) successions (e.g. Dominant → Tonic or Tonic → Subdominant).

The mechanisms are implemented in an independent way (when possible) and weighted using a standard range [0; 1], in order to easily configure (or eliminate) their degree of influence.

Data representation of harmonic functions

Each harmonic function is stored in a structure that contains a vector of Harmonic Function Strengths (since we consider 7 different diatonic functions its length is 7) and a corresponding vector of function modifiers. Harmonic Function Strength is a factor used to determine probability of occurrence of a specific diatonic harmonic function. In the authors' experiments the function modifiers are mostly limited to sevenths, with occasional ninths and in rare cases to sixths. Sixth is present in so-called Chopin's chord : Dominant with seventh and natural or augmented sixth instead of fifth, resolving down to the first degree of scale. The diatonic function with modifiers in a given key is equivalent to chords e.g. Tonic with modifier seventh in key G-major is equivalent to G 7 . Some modifiers cover others, e.g. a chord with a ninth is a variation of a chord with a seventh. Modifiers can be seen as subclasses of the main class (basic diatonic function).

For simplification reasons the authors have decided to determine the harmonic functions in constant intervals: twice per measure in the case of 3 or 4 beats and once per measure in the case of 2 beats. This is sufficient in most cases, as in practice the harmonic functions rarely change more frequently. The authors have also experimented with harmonic fragment being equal to the whole measure. Musical piece is decomposed into harmonic fragments (defined as indivisible musical unit with a single diatonic function attached to each one of them). The Harmonic Function Strength for a harmonic fragment is a sum of function excitations (determined by note pitch) for all notes in the fragment, according to varied degrees of Note Importance.

The flow of the proposed model

The flow of the proposed model regarding a single harmonic fragment is presented in Fig. 1, and proceeds as follows:

1. The harmony of a particular fragment is evaluated by examining all notes that it is composed of. Based on a set of rules from Section 3.3, every note is attached with a Note Importance value. 2. Each note excite (add a certain value to Function Strengths vector) several corresponding diatonic functions by being a specific chord component, with degree defined by the Function Excitation Matrix (Section 3.4). 3. After all the notes are processed, the vector of Function Strengths is elementwise multiplied by a Function Popularity Vector (Section 3.5). This primarily serves as a method of favouring the more popular functions. The Function Popularity Vector can also be a way to get a rich and uncommon harmony. 4. A Function Successions mechanism is applied after the individual processing of all of the harmonic fragments is complete. Function Succession matrix (Section 3.6) defines the degrees of desire for a specific succession of harmonic functions. They serve as additional modifiers (regarding direct predecessor and successor) for vectors of the Function Strengths. Every fragment is modified twice, first as a predecessor and the second time as a successor; only Each omittable mechanism (Note Importance, Function Popularity, Function Succession) is attached with an overall Global Configuration Weight from the range [0; 1] (Fig. 1). It can be used to fine tune or disable the influence of a specific mechanism on the final harmonization result. Excitation is an essential mechanism to obtaining Function Strengths, therefore it may not be weighted in a similar manner. The mechanisms are based on matrices (Function Excitation, Function Succession), vectors (Function Popularity) or rules (Note Importance) that offer substantial, direct and indirect configuration possibilities. The following subsections describe in detail the above-mentioned mechanisms.

Note importance determination

Note Importance may be determined by:

-the length of note -longer notes have, in general, greater influence on the harmony. The note length is used as an initial value for further importance weighting. Values corresponding to neighboring note lengths (alike eighth note and quarter note) have common ratio, empirically stated as 5/6. -position in measure -notes at inherently accented parts of measure -onbeat -are very important. Notes at generally unaccentuated parts of measure -off-beat, are less significant. Notes occurring in-between these main beats are harmonically even less important. -notes that are at the end of ties have minor contribution to the harmony; when played on string instruments (e.g. grand piano, upright piano, guitar) they are not hit again. They are therefore relatively quiet. -notes that are easily-heard by human ear are placed in extreme voices (highest and lowest notes); notes in the middle voices are slightly harder to hear and therefore contribute less to the local harmony; -accentuation increases the volume and therefore harmonic value of the note.

Homophony harmonization was chosen for the experimental evaluation, therefore the authors have determined a set of Note Importance rules described in Table 1. In this case, rules regarding voice position are irrelevant, and are not taken into consideration. These rules would be however important for music works with many voices. Corresponding note length from upper part of the Table 1 is taken as an initial value of Note Importance. Initial values for non-standard note lengths (e.g. a half note with a dot) are determined proportionally. The weight corresponding to the first fitting rule from lower part of the Table 1 serves as the multiplier for the initial value.

Implementation of Functions Excitation by notes

A well-defined diatonic function contains at minimum three components (typically root note, third and fifth). More complicated diatonic functions contain also additional or altered intervals. Determination of exact pitches of these function components requires detection of the piece key and reading of the current key signature.

Two matrices of Function Excitation are used to define the excitations of diatonic functions with notes (from the melody or the accompaniment). The matrix for major scales is presented in Table 2. Each diatonic function is excited by the particular matching note to a specified degree. This, obviously, is relative to the pitch of the note in question, and to the pitches of the notes occurring in the tonal functions themselves. The authors have defined the weights of excitation in the range of [0, 1], assigning:

-[0.8; 1.0] -greater degrees to the common chord components (e.g. root note, third and fifth); -[0.4; 0.6] -moderate degrees to the less common components: sixth, seventh, ninth with occasional modification of natural third. Making major chords from naturally minor chords using augmented third serves frequently as tonicization -a local Dominant → Tonic resolution); -[0.1; 0.2] -small degrees to the rare modifications of natural sixth seventh and ninth. 0 0 0.9 0 0 0.8 0 0 0 0 0 T6 1 0 0 0.9 0.1 0.1 0 0.8 0 0 0.4 0 D 1 0.5 0.2 0.1 0.9 0.1 0 0.8 0 0.4 0.6 0.1 S 1 0 0 0.1 0.9 0.1 0 0.8 0 0 0.4 0 T3 1 0 0 0.9 0.1 0.1 0 0.8 0 0 0.4 0 S2 1 0 0 0.9 0.1 0.1 0 0.8 0 0 0.4 0 T 1 0 0 0.1 0.9 0.1 0 0.8 0 0 0.4 0

A similar matrix has been prepared for minor scales. Due to space limitation it is not presented in this paper. The obvious differences in reference to major scales are the natural qualities of diatonic functions. The less straightforward dissimilarity is a frequent conversion from naturally minor Subdominant and Dominant to major. This commonly occurs in a melodic variation of scale, and sometimes, in harmonic variations.

Function popularity

The goal of determining and applying weights corresponding to Function Popularity is to directly prioritize frequently occurring tonal functions (e.g. Tonic, Subdominant, Dominant). This is a direct way to make them more frequent than less popular functions, as inherently specified by music theory. The goal can be indirectly obtainable by using lower coefficients in the Functions Succession matrix (described in Section 3.6). Obviously, the direct method makes controlling the process much easier. For experimental studies the authors have used the Function Popularity vector, specified in Table 3.

Succession of harmonic functions

The succession of harmonic functions (horizontal relations between neighboring harmonic functions) is implemented using encouragement of more probable combinations, with moderation of less likable ones. This is done using a Function Succession matrix. The matrices proposed for major and minor scales are presented in Table 4. The matrix for minor scale is similar and due to space limitations not presented here. The matrices values were determined in order to prioritize the most likable successions. The exemplary common successions are: Tonic into Subdominant, Subdominant into Dominant and Dominant into Tonic -cadence (the simplest and most common). Another example is less common deceptive cadence (Dominant into T onic 6 ). The matrices also support tonicization (e.g. succession from Subdominant 2 into Dominant). It is worth mentioning that the most supported

Table 4. Functions Succession weights for major scales

Transfer into: T S2 T3 S D T6 D7 D7 0.5 0.3 0.5 0.2 0.9 0.1 1.0 T6 0.4 0.8 0.3 0.8 0.7 1.0 0.2 D 0.9 0.4 0.7 0.3 1.0 0.8 0.3 S 0.7 0.4 0.4 1.0 0.9 0.5 0.1 T3 0.3 0.2 1.0 0.5 0.6 0.8 0.1 S2 0.2 1.0 0.2 0.5 0.9 0.4 0.2 T 1.0 0.4 0.3 0.8 0.9 0.6 0.1 quasi-succession is maintaining the current function (no change) with the maximum degree of support: 1.0. It allows more efficient handling of common cases, where harmonic functions change less frequently than the arbitrary harmonic fragment length (e.g. two beats). It is also important to mention that the support of succession occurs only between harmonically determinable harmonic fragments, excluding these that do not contain notes at all (only pauses), or contain only ongoing tied notes from previous beats. It is assumed that such fragments continue the previous harmonic function (which is a slight oversimplification as it does not have to be always true). The succession support occurs forwards and backwards for each possible succession, with a degree defined by the sum of products of the neighboring Function Strengths and a succession weight from the Functions Succession matrix.

Model tuning and experimental results

The authors have implemented the proposed model and tested it using music documents in MusicXML file format [START_REF] Good | MusicXML for Notation and Analysis[END_REF]. This section discuses various tuning possibilities and presents exemplary homophony harmonizations with the use of various tuning approaches.

The model has been applied to musical works containing a single melody (monophony) with no key changes (no modulations). It may be considered as a greater challenge (more harmonic uncertainty) than using musical works with accompaniment or several independent voices, where harmonic functions are generally easier to detect. The proposed harmonization model is viable for almost every musical piece, regardless of the number of voices, as long as it is maintained in a tonal system, e.g. uses either major or minor scale (with possible modifications like harmonic, melodic, Dorian, etc.). In case of many voices or chords, ideally, a customized Notes Importance determination is required (as described in Section 3.3), in order to detect and prioritize more important voices, and attenuate the less important ones. Efficient harmonization of a musical works with key changes (modulations) requires detection or indication of such changes, and an adequate update of tonal root and/or minor/major scale properties.

The authors propose the following model modifications: The learning set has relatively simple harmonies, therefore the harmonization results (produced by the tuned model) are also simplified.

Regarding indirect tuning: it is possible to tune the original model into a specific direction by training with a preselected set of examples (i.e. jazz pieces with complicated harmony or musical pieces from a specific period). Such approach requires a number of carefully selected examples and is problematic as musical styles are often ambiguous.

Remarks on universality

The proposed model offers numerous configuration possibilities and universality: as it is not data-driven but knowledge-driven, it provides multiple levels of model control and, to some extent, results of harmonization. As opposed to the data-driven approaches (taught with examples) the authors' methodology does not need to be fed with large or representative sample of data (requiring data gathering, selection and frequently producing uncontrollable results). It may be used in many variants relying only on the theoretical knowledge. With modifications of underlaying configuration data one can achieve many valuable tasks, such as:

defining the Note Importance rules tuned to various musical pieces; defining the Function Excitation matrix to generate rare or simple functions/modifications, resulting in a complicated or simplistic harmony; defining the Function Popularity vector to moderate, forbid or encourage specific tonal functions, directly limiting or increasing probability of their occurrence; defining the Function Succession matrix to moderate, forbid or encourage specific harmonic function successions; changing the Global Configuration Weights of harmonization mechanisms in order to tune the model to the specific needs and expectations; iterate through configurations in order to quickly generate various (possibly interesting) harmonies for the same musical work.

The authors propose also data-driven tuning of the initial model, what could provide more viable harmonizations or (with preselected set of examples) produce specific harmony style.

Future works

Future works in the area will be conducted in the following directions:

continuous development and evaluation of the presented harmonization model; determination of configuration matrices and parameters for different styles of musical works (like jazz, classical music, popular music, etc.) using indirect tuning with preselected set of examples (what requires relatively large and representative examples database); development of the Function Succession mechanism with use of chord modifiers; development of the Function Excitation mechanism using also relations between notes, rather than independent notes; developing criteria for automatic evaluation of the obtained harmony; automatic parametrization of the harmonization model based on the abovementioned criteria; further development of model tuning procedures, using well established machine learning paradigms, e.g. Artificial Neural Networks.

Fig. 1 .

 1 Fig. 1. Harmonization model flow

Fig. 2 .

 2 Fig. 2. Alternative harmonizations using different tuning procedures

Table 1 .

 1 Determining Note Importance

	Note length	Initial value
	Whole note	1.440
	Half note	1.200
	Quarter note	1.000
	Eighth note	0.833
	Sixteenth note	0.694
	Thirteenth-second note	0.579
	Condition	Importance weight
	Note at the end of a tie	0.4
	Note on the 1 st beat	1.0
	Note on the 3 rd beat	0.9
	Accentuated note	0.8
	Note on the 2 nd or 4 th beat	0.7
	Other notes	0.6

Table 2 .

 2 Excitation weights for major scales

	Absolute shift from function root note[in semitones]
	0	1	2	3	4	5	6	7	8	9 10 11
	D7 1

Table 3 .

 3 Function Popularity vector weights

	Tonal function
	Name	Weight
	T	1.0
	S2	0.6
	T3	0.7
	S	0.9
	D	0.95
	T6	0.65
	D7	0.2

  Function Popularity vector and Function Succession matrix) using non-direct modification (tuning) of matrices. The example is a harmonized music piece: a music document with attached harmony functions corresponding to class labels. A simple statistical analysis of the occurring harmonic functions are used to update the values in the Function Popularity vector. Analysis of occurring harmonic function successions in the example are used to update Function Succession matrix. Total extinction of rare functions and successions may occur (depending on the harmony relations included in the examples), therefore the authors propose to limit the minimal values to a fixed threshold 0.05. The learning process is independent regarding the sequence of examples -similarly to batch learning, as sequential learning would prioritize the recent examples over the previous ones. {0, 0, 0.8} in the second case. In the second case less popular functions are encouraged to appear. Third harmony was obtained by using direct tuning of the model (by changing Function Excitation matrix to produce more seventh and ninth chords, and Function Popularity vector to encourage rare functions. Fourth harmony was obtained by using indirect tuning by learning from examples, with ten examples being modern popular music pieces of unspecified genre.

	1. changing values of Global Configuration Weights (respectively: Note
	Importance, Function Popularity, Function Succession), exemplary applica-
	tions are: moderation of less popular functions, introduction of uncommon
	function successions;
	2. direct tuning -changing values in the configuration matrices (Function
	Excitation matrix, Function Popularity vector, Function Succession matrix)
	and defining custom Note Importance determination rules. The initial values
	are meant to be universal, changing them directly influences the behavior of
	the model and produced results in a specific direction, but requires experi-
	ence and advanced knowledge of harmony;
	3. indirect tuning by learning from examples -modification of a part
	of proposed approach (The results in Fig. 2 present alternative harmonizations produced using different
	tuning procedures. The fragment of musical piece (A. De Vita and H. Sharper -
	Softly, as I leave you) was harmonized with harmonic fragments set to 2 beats
	(half of measure). The first and second harmonies were obtained for original
	model with various values of Global Configuration Weights (respectively: Note
	Importance, Function Popularity, Function Succession): {0.4, 0.4, 0.4} in the first
	case,

Acknowledgement

This work is supported by The National Center for Research and Development, Grant no. N R02 0019 06/2009.