Proceedings of the 2014 eC 1st Student Computer Science Research oSR Conference C Stu University of Primorska Press StuCoSReC Preface Proceedings of the 2014 1st Student Computer Science Research Conference Computer science is one of the quickest developing ar- eas of Science. This area captures a lot of disciplines Edited by that have been developed independently. Thus, new dis- Iztok Fister jr. and Andrej Brodnik ciplines have been arisen. On the other hand, the com- puter science starts connecting with the other natural Reviewers and Programme Committee sciences in order to draw an inspiration for solving their Janez Brest, Chair ■ University of Maribor, Slovenia Andrej Brodnik, Chair problems. Nowadays, incorporating principles from biol- ■ University of Primorska and University of Ljubljana, Slovenia ogy (e.g., Darwinian evolution or behavior of social liv- Iztok Fister, Chair ■ University of Maribor, Slovenia ing insects and animals) into the computer algorithms Iztok Fister jr., Chair ■ University of Maribor, Slovenia have been observed. In general, this means enough ma- Zoran Bosnić ■ University of Ljubljana, Slovenia terial for a computer science conference. This volume Borko Bošković ■ University of Maribor contains papers as presented in the first Student Com- Mojca Ciglarič ■ University of Ljubljana, Slovenia puter Science Research Conference 2014 ( StuCoSRec) Simon Fong ■ University of Macau, Macau held in Maribor under the cover of the 17th MultiCon- Tao Gao ■ North China Electric Power University, China ference on Information Society in Ljubljana and ACM Marjan Heričko ■ University of Maribor, Slovenia Slovenia. Andres Iglesias ■ Universidad de Cantabria, Spain Branko Kavšek ■ University of Primorska, Slovenia The 1st Student Computer Science Research Conference Miklos Krecz ■ University of Szeged is an answer to the fact that modern PhD and MSc pro- Gregor Papa ■ Jožef Stefan Institute grams foster early research activity among the students Peter Rogelj ■ University of Primorska, Slovenia of computer science. The prime goal of the conference is Xin-She Yang ■ Middlesex University, United Kingdom to become a place for students to present their research Borut Žalik ■ University of Maribor, Slovenia work and hence further encourage students for an early Organizing Committee research. Besides the conference, it also wants to estab- Uroš Mlakar ■ University of Maribor, Slovenia lish an envi ronment where students from different insti- Jani Dogonik ■ University of Maribor, Slovenia tutions meet, let know each other, exchange the ideas, Published by and nonetheless make friends and research colleagues. University of Primorska Press In line with this, we would like to exhibit research ef- Titov trg 4, si-6000 Koper forts of students in four Slovenian institutions, i.e., Uni- Editor-in-Chief versity of Maribor, University of Ljubljana, University Jonatan Vinkler of Primorska and Jožef Stefan international post-gradu- ate school, neighboring countries, i.e., Austria, Croatia, Managing Editor Hungary and Italy, and the other countries in the world Alen Ježovnik (e.g., one of the papers on out conference comes from Koper, 2014 China). At last but not least, the conference is also meant to serve as meeting place for students with senior isbn 978-961-6963-03-9 (pdf) www.hippocampus.si/ researchers from different institutions. isbn/978-961-6963-03-9.pdf isbn 978-961-6963-04-6 (html) Eleven papers addressed this conference,1 covering sev- www.hippocampus.si/isbn/978-961-6963-04-6/index.html eral topics of the computer science. All the papers were © 2014 Založba Univerze na Primorskem reviewed by two international reviewers and accept- ed for the oral presentation. This fact confirms a good work with authors in their research institutions. The content of the papers will be presented in three sections that cover theory and applications, computer graphics, image processing and pattern recognition, and semantic web and data mining. The conference is dedicated to graduate and un- der-graduate students of computer science and is there- CIP - Kataložni zapis o publikaciji fore free of charge. In line with this, the organizing Narodna in univerzitetna knjižnica, Ljubljana committee would like to thank to the University of Maribor, especially the Faculty of Electrical Engineer- 004(082)(0.034.2) ing and Computer Science (FERI) for the support. Es- STUCOSREC [Elektronski vir] : proceedings of the 2014 pecially, we are grateful to the Institute of comput- 1st Student Computer Science Research Conference / edited er science at FERI Maribor for payment of conference by Iztok Fister jr. and Andrej Brodnik. - El. knjiga. - Koper : costs arising during the organization. At the end, a spe-University of Primorska Press, 2014 cial thank goes to the dean of UM FERI, prof. Borut Način dostopa (URL): www.hippocampus.si/isbn/978-961-6963- Žalik for his unselfish support. 03-9.pdf Način dostopa (URL): www.hippocampus.si/isbn/978-961-6963- 04-6/index.html ISBN 978-961-6963-03-9 (pdf) ISBN 978-961-6963-04-6 (html) 1. Fister, Iztok, ml. 275792896 1 Conference site is at http://labraj.uni-mb.si/stucosrec2014/. Contents Preface II Spatially Embedded Complex Network Estimation Using Fractal Dimension ■ David Jesenko, Domen Mongus, Borut Žalik 5–8 Graph Coloring Based Heuristic for Driver Rostering ■ László Hajdu, Miklós Krész, Attila Tóth 9–12 A New Heuristic Approach for the Vehicle and Driver Scheduling Problems ■ Viktor Árgilán, Balázs Dávid 13–17 Multi Population Firey Algorithm ■ Jani Dugonik, Iztok Fister 19–23 Monte Carlo Path Tracing with OpenCL ■ Andrej Bukošek 25–27 A Garlic Clove Direction Detection Based on Pixel Counting ■ Pavel Fičur 29–31 Simple Approach to Find Repeated Patterns in Opponents Texas Hold‘em No-limit Game ■ Gregor Vohl, Janez Brest 33–36 Histogram of Oriented Gradients Parameter Optimization for Facial Expression Recognition ■ Uroš Mlakar 37–41 Constructing Domain-specific Semantic Dictionaries to Supplement Domain-specific Knowledge Bases ■ Goran Hrovat, Milan Ojsteršek 43–45 Am I Overtraining? A Novel Data Mining Approach for Avoiding Overtraining ■ Iztok Fister Jr., Goran Hrovat, Samo Rauter, Iztok Fister 47–52 An Improved Algorithm of DPM for Two-dimensional Barcode ■ Tao Gao, Xiao-cheng Du, Iztok Fister Jr. 53–56 StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October III StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October IV Spatially Embedded Complex Network Estimation Using Fractal Dimension ∗ David Jesenko Domen Mongus Borut Žalik University of Maribor University of Maribor University of Maribor Faculty of Electrical Faculty of Electrical Faculty of Electrical Engineering and Computer Engineering and Computer Engineering and Computer Science Science Science Smetanova 17, SI-2000 Smetanova 17, SI-2000 Smetanova 17, SI-2000 Maribor, Slovenia Maribor, Slovenia Maribor, Slovenia david.jesenko1@um.si domen.mongus@um.si borut.zalik@um.si ABSTRACT Peitgen et al. introduced three main Mandelbrot fractal This paper describes complex networks in regards to their dimensions [4]: fractal characteristic. Firstly, a short background concern- ing fractal dimensionality based on box-counting is consid- ered, while in the continuation its application on spatially • self-similarity dimension, embedded networks is proposed using rasterization at mul- • compass dimension (also called divider dimension), tiple scales. Finally, the obtained results are presented. • box-counting dimension. General Terms Computer science, algorithms The box-counting dimension is the most frequently used in computer-supported applications. The box-counting algo- Keywords rithm uses a uniform grid with cell-size s1, where the num- Complex networks, fractals, fractal dimension ber of grid-cells N1 containing the object is counted. Obvi- ously, at an increased sampling scale s2, an increased value 1. INTRODUCTION of grid-cells containing structures N2 is expected. A fractal dimension D, is then given by [5]: Today, complex networks are amongst the most advanced approaches for describing and analyzing structures composed log N1 − log N2 D = . (1) of a plurality of blocks. When vertex has a fixed position in log s1 − log s2 space, we are talking about spatially embedded networks [1]. Measuring fractal dimension with the box-counting algo- Similar to graphs, networks are represented by a set of ver- rithm has been presented at various places. Conci and Proen- tices and an unordered set of edges linking them. The main ça in [6] used this method for visual inspection, which is difference between them is of topological nature. Complex important part of quality control in textile industry. In [7] networks have nontrivial and unpredictable topology [2]. In Boshoff presented fast box-counting algorithm for determin- addition, topology of complex networks hide important in- ing the fractal dimension of sampled continuous functions. formation about vertices. Evaluation of these features can Zhang et al. in [8] presented coarse iris classification using reveal many otherwise hidden characteristics of geometric box-counting algorithm to estimate fractal dimensions. structures as for example separation, connectedness, com- pactness, or countability conditions [3]. The actual fractals are difficult to be detected and evaluated. Frequently, a definition of fractal dimension is used for this With their characteristic self-similar property and huge amo-purposes [4]. The dimension of fractals is usually smaller unt of possibilities to analyse them, fractals are one of the than the space in which fractals are embedded. However, most interesting areas of computer science. Among various the fractal dimension remains the same regardless of the res- dimensions that are extensively studied in the literature, olution used for fractal dimensionality estimation. Although ∗Corresponding author. this definition allows for identification of fractal geometry, fractal dimension of topology is still not clearly defined [9]. In this paper we present an implementation of a box-counting method for estimating fractal dimension of topology of spa- tially embedded complex networks. Section 2 gives a short theoretical background on complex networks. In Section 3, a method for estimating fractal dimensions of topologies is presented. Implementation details are given in Section 4. Results are presented in Section 5, while Section 6 concludes the paper. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 5 2. COMPLEX NETWORKS Nature consists of many complex networks established over various kinds of vertices with different relations. Obvious examples are social networks, where vertices represent peo- ple and edges defines different relations between them, e.g. business, friendship, or family [10]. In all these cases, a network N is defined by the set of vertices V (N ) = {vi} and the Figure 3: Koch snowflake. set of edges E(N ) = {eij }, where eij = (vi, vj ) represents a pair of vertices vi, vj ∈ V (N ). 2, two types of spatially embedded networks exist. An ex- 2.1 Spatially embedded complex networks ample of a network where vertices are connected with their Depending on the length of the edges between vertices, two nearest neighbours is shown in Figure 4. For creating these kinds of spatially embedded networks are known [11]: edges types of networks, an efficient method for searching near- can either exist only between the nearest neighbours, or they est neighbours, implemented in nano-flann library, available may spend over the entire network. In all the cases, power in [12] has been used. It is based on approximate nearest law distribution of edge lengths is given as [11]: neighbours search algorithm introduced in [13]. In Figure 4, the brighter vertices are those with more edges than the darker one. P (r) ∼ r−δ, (2) where r is the distance between the vertices. Exponent δ affects the network dimension as: • if 4 < δ, dimension is equal 2, because the value of P (r) is low. • if 2 ≤ δ ≤ 4, dimension monotonous grows from 2 → ∞ with decreasing δ. • if δ < 2, dimension converts towards infinity. 3. FRACTALS Fractals are geometric objects, which cannot be described Figure 4: Complex network containing edges to by ideal or primitive geometric objects, like cubes, spheres, their nearest neighbours. triangles, or polygons. Usually they are related with recur- sive operations on geometric objects and number sets [4]. When considering complex networks (as shown in Figure 5) These operations are mathematically straightforward, how- with edges spanning over the entire space, the implemen- ever the results are complex, with interesting and remark- tation was based on similarity measurement of vertices (ac- able properties. The most important property of fractals is cording to a user-specified feature). Obviously, the construc- self-similarity [4]. Examples of this can be found in the Can- tion of this type of network is significantly more expensive as tor set (Figure 1), Sierpinski triangle (Figure 2) and Koch it demands comparison of each pair of vertices. This leads to snowflake (Figure 3). a quadratic time complexity, while approximate search for the nearest neighbours can be achieved in logarithmic time [14]. Figure 1: Cantor set. Figure 2: Sierpinski triangle. 4. IMPLEMENTATION First, an application for generating different kind of com- Figure 5: Complex network with edges over the plex networks has been developed. As explained in Section whole network. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 6 4.1 Box-counting algorithm A double-step algorithm using only integer arithmetic was used for the implementation of line rasterization with min- imized oversampling. As proposed in [15], the method uses the coding of changes in directions based on DDA (Digi- tal Differential Analysis). In order to increase efficiency of box-counting, the scales of both uniform superimposed grids should be in the ration 1:2. In this approach each box from a grid is subdivided into four boxes each of half the size in the next grid. Equation for fractal dimension is now [4]: log N (2−(k+1)) − log N (2−k) N (2−(k+1)) = log . (3) log 2k+1 − log 2k 2 N (2−k) Figure 7: Fractal superimposed with grid, where s = 1/40 width of superimposed grid. The result uses the base 2 logarithm due to the factor by which the box-count increases from one grid-level to the next [4]. of |V (M )| = 250 vertices and |E(M )| = 2554 edges. With an increasing scale, CPU time is decreased. As seen, the box- 5. RESULTS counting fractal dimension is practically stable. The varia- Implemented algorithm runs on a computer system with the tions can be assigned to the point-to-cell algorithm, which following configuration: approximate the complex network more accurately at larger resolutions. • processor Intel Core2 Quad CPU Q9550 with beat 2.83GHz, Table 1: Box-counting fractal dimension. • memory DDR3 of size 8GB, Cell size CPU time [ms] Fractal dimension 1/35 and 1/40 59.4 1.31 • graphical card ATI Radeon HD 4600 Series with 1GB 1/30 and 1/40 56.1 1.33 video memory, 1/25 and 1/40 50.4 1.33 • 1/20 and 1/40 45.3 1.34 operating system Windows 7 Pro - 64 bit. 1/20 and 1/35 39.7 1.35 1/20 and 1/30 33.5 1.35 Figures 6 and 7 shows spatially-embedded complex networks 1/20 and 1/25 28.9 1.36 superimposed on grids of two different scales. For practical 1/15 and 1/20 26.2 1.37 purposes it is often convenient to consider a sequence of grids where the mesh size is reduced by a factor of 1/2 from one grid to the next. In other words, if the number of boxes The spent CPU efficiency of the implementation regarding counted increases by a factor of 2D when the box size is different network size is given in Tables 2 and 3. The scales halved, then the fractal dimension is equal to D. More in used were 1/20 and 1/40 width of superimposed grid. As [4]. expected, the larger networks demand longer time. Differ- ences in CPU efficiency between both kinds of networks can be noticed when comparing Tables 2 and 3. It is a conse- quence of rasterization algorithm. Edges in networks, with edges over the whole network, are in average longer and their rasterization therefore take more time. Table 2: Spent CPU time for evaluation of spatial dimension with nearest neighbour edges. Number of vertices Number of edges [millions] Time [s] 100 0.001 0.01 500 0.008 0.1 1000 0.04 0.3 5000 0.9 0.6 10000 3.6 1.1 Figure 6: Fractal superimposed with grid, where s = 25000 23.7 4.7 1/20 width of superimposed grid. 50000 92.1 18.2 75000 207.5 37.4 In Table 1, fractal dimension calculated by Eqn. (1) at a changing scale can be seen. Complex network was composed StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 7 https://code.google.com/p/nanoflann. Table 3: Spent CPU time for evaluation of spatial [13] S. Arya, D. M. Mount, N. S. Netanyahu, R. Silverman, dimension with edges over the whole network. and A. Y. Wu. An optimal algorithm for approximate nearest neighbor searching fixed dimensions. Journal Number of vertices Number of edges [millions] Time [s] of the ACM (JACM), 45(6):891–923, 1998. 100 0.0008 0.015 500 0.005 0.2 [14] J. L. Bentley. Multidimensional binary search trees 1000 0.03 0.45 used for associative searching. Communications of the 5000 0.8 0.7 ACM, 18(9):509–517, 1975. 10000 3.1 1.3 [15] Yong Kui Liu, Peng Jie Wang, Dan Dan Zhao, Denis 25000 18.3 5.4 Spelic, Domen Mongus, and Borut Zalik. Pixel-level 50000 79.8 23.8 algorithms for drawing curves. In Theory and Practice 75000 176.1 47.1 of Computer Graphics, pages 33–40. The Eurographics Association, 2011. 6. CONCLUSION In this paper, a brief theoretical background on complex networks and fractals (fractal dimension) has been given at first. The implementation of the box-counting algorithm for evaluation of fractal dimension of spatially embedded com- plex networks follows. The efficiency of the implementation is given for two types of complex networks. 7. REFERENCES [1] L. Barnett, E. Di Paolo, and S. Bullock. Spatially embedded random networks. Physical Review E, 76(5):056115, 2007. [2] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. U. Hwang. Complex networks: Structure and dynamics. Physics reports, 424(4):175–308, 2006. [3] R. Engelking. General topology. 1989. [4] H.O. Peitgen, H. Jürgens, and D. Saupe. Chaos and fractals: new frontiers of science. Springer, 2004. [5] B.S. Raghavendra and N.D. Dutt. Computing fractal dimension of signals using multiresolution box-counting method. International Journal of Information and Mathematical Sciences, 6(1):50–65, 2010. [6] A. Conci and C. B. Proença. A fractal image analysis system for fabric inspection based on a box-counting method. Computer Networks and ISDN Systems, 30(20):1887–1895, 1998. [7] H. F. V. Boshoff. A fast box counting algorithm for determining the fractal dimension of sampled continuous functions. In Communications and Signal Processing, 1992. COMSIG’92., Proceedings of the 1992 South African Symposium on, pages 43–48. IEEE, 1992. [8] L. Yu, D. Zhang, K. Wang, and W. Yang. Coarse iris classification using box-counting to estimate fractal dimensions. Pattern Recognition, 38(11):1791–1798, 2005. [9] O. Shanker. Defining dimension of a complex network. Modern Physics Letters B, 21(06):321–326, 2007. [10] M. Aldana. Complex networks. Available on: http://www.fis.unam.mx/ max/English/notasredes.pdf, 2006. [11] L. Daqing, K. Kosmidis, A. Bunde, and S. Havlin. Dimension of spatially embedded networks. Nature Physics, 7(6):481–484, 2011. [12] Nano-flann. Available on: StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 8 Graph coloring based heuristic for driver rostering László Hajdu Miklós Krész Attila Tóth University of Szeged University of Szeged University of Szeged Dugonics square 13. Dugonics square 13. Dugonics square 13. Szeged, Hungary Szeged, Hungary Szeged, Hungary hajdul@jgypk.u- kresz@jgypk.u- attila@jgypk.u-szeged.hu szeged.hu szeged.hu ABSTRACT means that the employees don’t spend the same amount of Nowadays, the companies and institutions have numerous time working. These differences may produce an overtime employees therefore the crew rostering problem became in- cost which is added to the basic salary. However, the com- creasingly important. The objective is to assign the crew panies and institutions must guarantee the basic salary for members to previously generated daily shifts, while meeting everyone, even if the employee does not spend the normal the constraints, and to optimize the overall cost. In this pa- amount of time at work. This causes an additional cost for per we present the mathematical model of the problem and the companies. introduce a two-phase graph coloring method for the crew rostering problem. Our method has been tested on artifi- In this paper we define our crew rostering problem, intro- cially generated and real life input data. The results of the duce the mathematical model and give our new heuristic new algorithm have been compared to the solutions of the solution for the problem. The results of the algorithm have appropriate integer programming model for problems with been compared to the results of the appropriate integer pro- moderate size. gramming model for problems with moderate size. These results show that our algorithm is very efficient for solving Categories and Subject Descriptors this problem. Our solution is provided for the driver ros- H.4 [Information Systems Applications]: Miscellaneous; tering problem, although it can be utilized in crew rostering H.4.2 [Types of Systems]: Decision support(e.g., MIS) within a wider range of propositions. General Terms Algorithms, Management, Measurement, Performance 2. CREW ROSTERING The problem is based on generalized set covering model. Keywords Dantzig was the first to deal with its mathematical appli- Driver rostering, Graph coloring, Tabu Search cation[4]. The literature gives numerous examples of the issue. The most significant ones are the airline [5], the call 1. INTRODUCTION center[12], the nurse scheduling [3] and the driver schedul- ing, what is the topic of our research. Since giving an ex- The crew rostering appears in various areas. The main chal- act solution to these problems is only possible by moderate lenge is to assign an optimal number of crew members to the sample size, thus these are solved usually with heuristics. pre-defined daily shifts in a way where the resources are used Several solutions can be found in the literature, such as: ant in the most efficiently way, meanwhile, the solution meets colony optimization, dynamic programming, genetic algo- the defined constraints of the workplace. It can be stated rithm, simulated annealing or tabu search method. A good that it is possible to use the model in every sector which overview of the above methods with respect to driver ros- deals with crew members. Some of the typical application tering can be found in [11]. areas where the crew rostering can be applied: public trans- port companies, call centers, and nurse scheduling. The vari- Our crew rostering problem is formally defined as the follow- ations of the problem are NP-hard or NP-complete therefore ing. Let c be the set of employees. A set of shifts denoted extremely hard to solve. [8, 1, 10] by S needs to be carried out. A shift is composed of a series of tasks. A shift is defined by its date, starting and ending In most cases the length of the shifts are different, which time, duty time, and its working time. The working time can be different from the duty time because of breaks or idle ac- tivities. Our objective is to assign the crew members to the shifts. Consequently let f = s → c an assignment where the shifts are covered by the employees in a way where there is a person assigned to every shift. These solutions need to correspond the regulations of the company as well as the norms of the european union. For example, that the crew members can get at most one shift in one day. Besides this, the following attributes were regulated during this problem: StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 9 • Minimum rest time between two shifts. Driver i gets shift k on day j: • Maximum worktime during one week. xijk ∈ {0, 1} • Minimum free days in one month. Driver i works on day j: • Maximum consecutive workdays. zij ∈ {0, 1} Driver i works in the planning period: The aim is to give a solution which meets the regulations yi ∈ {0, 1} above, while minimizing the cost. Let ct(i) be the type of contract belonging to driver i. For any contract type c, let Overtime of driver i: aw(c) return the expected daily working hours that belong π to c. Each driver i has such a contract that defines his i ≥ 0 required daily working time in average. Based on their con- tract it is possible to calculate the expected worktime in the The constraints of the appropriate integer programming model planning period. In our case every crew member had a con- were the followings: tract with eight hours of worktime. The expected worktime could be defined the following way: Every driver needs to be assigned to exactly one shift. ∀j∀k xijk = 1 (1) i expected worktime = number of workdays ∗ aw(ct(i)) A driver only works on a certain day if there is a shift as- signed to them and only one person can be assigned to a The employee needs to be paid based on their expected work- shift. time written in their contract. Also, in the case of working over the expected worktime, the employers have to pay extra ∀i∀j xijk = zij (2) for this overtime. The optimal case is when every employee i works in their expected worktime. Therefore, the cost is the A driver is employed if they are working in the planning following: period. ∀i zij ≤ Lyi (3) j cost = α ∗ overtime + β ∗ employment cost The following constraint excludes the possibility of assigning two incompatible shifts to a driver. Where α, β are weights. Hence the objective was that the overtime is as low as possible while the number of the em- ∀i∀j∀p∀k∀qxijk + xipq ≤ SSjk pq + 1 (4) ployees is minimal. The worktime should not exceed the maximum working time in a week. 3. MATHEMATICAL MODEL Let D be the days of the planning period, DW eek the days of ∀i∀w xijkwt(j, k) ≤ W T (5) the week and DMon is the days of a months. Meanwhile the W k j∈D eek W lenght of the planning period needs to be defined denoted by L, and let L A driver can only have the maximum consecutive shifts. m the number of days in the months denoted by m. Furthermore let j and p (j ∈ D, p ∈ D) the days ∀i∀p z and w the weeks. The set of the drivers are denoted by C, ij ≤ W D (6) and let S be the shifts where S j∈DW D j is a shift on day j. Let i p (i ∈ C) the driver i and k and q (k ∈ Sj , q ∈ Sj ) one shift At least the required free days need to be given in a month. each. SSjk pq is a compatibility relation, which takes up the value of 1 if the shift k on day j and shift q on day p can ∀i∀m zij ≤ Lm − RD (7) be assigned to the same person, and 0 otherwise. Let W T j∈DM on be the maximal worktime on a week, W D the maximum m number of consecutive workdays and DW D Let define the overtime. p the consecutive days beginning from day p (W D+1). The maximum number ∀i ( x of free days is denoted by RD. For every driver a worktime ijk wt(j, k) − ct(i, m)) ≤ πi (8) being defined in the contract is needed. Let this worktime m j∈DM on k m be ct(i, m). In order to be able to minimize the number of employees let define cc(i) as the operational cost. Let The objective function minimizes the sum of the overtime wt(j, k) be and st(j, k) be the worktime and the duty time and the operational cost. of every shift k on every day j. Finally let ET be the value used by multiplying the overtime. We needed the followings in order to be able to define the min (ET πi + yicc(i)) (9) model: i StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 10 4. HEURISTIC METHOD the shifts every person’s worktime could be calculated for The algorithm is a two phase heuristic method based on the initial solution. This way, since the real working time graph coloring. M. Gamache et. al [7] developed a method, of each driver and their expected worktime were known, an where graph coloring algorithm based tabu search was used initial cost could be calculated for the solution. for scheduling pilots. This gave the idea for our problem where the rules being discussed above had to be met. Be- 4.3 Tabu Search sides the solution had to match the important aspects of the The tabu search was introduced by Glover in 1986 and it driver rostering. The heuristic has the following steps: is one of the most famous local search techniques nowadays [9]. The state space of the tabu search consists every feasi- 4.1 Initial steps ble coloring, and every state is a coloring. Let SL be all of 4.1.1 Set of employees the possible solutions of the graph coloring. The objective function had to be minimized in SL. Every sl ∈ SL pos- We can give an estimation for the number of employees for sible solution has a neighbourhood denoted by N (sl). The the problem using the worktime of the shifts and the con- algorithm visits sl tracts of the workers. Using this, a lower bound for the 0, sl1, ...sln solutions during the runtime where sl total overtime can be calculated. The lower bound of the 0 is the initial solution produced by the graph col- oring, and sl algorithm depends on the number of the employees and the i+1 ∈ N (sli). The next neighbour is chosen by first fit method. This means that the first neighbour solution working time. The number of the crew members is the fol-being better than the actual solution is chosen. Neighbour- lowing: hoods chosen in one step are stored in the tabu list denoted by T L. The steps in the list are not being performed for a time. The N (sl) set of neighbours of the sl solution can be global worktime C = round( ) defined by using the following neighbourhoods: contract type 1. Recoloring of a vertex: During this step a new color is It is important to note that a trivial lower bound for the given to a vertex of the graph. Meaning that a shift number of employees can also be given by the number of the is taken a way from a driver and given to another one shifts on the busiest days, since one crew member can have who is able to perform the task while keeping the con- at most one shift a day. straints. Shifts are taken away from the drivers having the most overtime and given to the ones having under- 4.1.2 Days-off patterns time. A days off patterns defines the fixed free days for an em- 2. Swapping the colors of two vertices: During the swap, ployee. For scheduling the free days of the employees, in the colors of two vertices are switched, therefore shifts many cases days off patterns are used for the heuristic solu- between two drivers are swapped in a way where both tions of crew rostering [6]. We propose a 6-1 days-off pattern, of them receives a shift which they can carry out keep-which means that the workers get one day off after every six ing the constraints. Usually this happens between the workdays. This pattern meets both of the minimal free days drivers having the most overtime and the ones having and the maximal consecutive workdays, and causes only a undertime. few exclusion during the graph coloring, and the tabu search will have a relatively big state space. After carrying out the vertex recoloring or the vertex swap- 4.1.3 Graph Building ping, the step is added to the tabu list. Multiple methods can be applied as a stopping criteria. In most cases the prac- Let every shift in the planning period be a vertex of a graph. tice uses criteria based on iteration number or time bound. There is an edge between two vertices if the shifts repre- In our case a stopping criteria was defined where the al- sented by the two linked vertices cannot be performed by gorithm stops if it can not find a better solution within a the same person. This way every day is a clicque and there certain number of iteration. are additional edges between the clicques because of the rest time constraint. 5. TEST RESULTS 4.2 Graph Coloring The solutions of the heuristic method were compared to the appropriate integer programming’s solutions using a rela- During graph coloring each color mean one of the drivers. tively small input (a maximum daily shifts of 50). The algo- The initial estimation was usually correct but some input rithm was implemented in Java language and IBM Cplex with too short of a worktime might have result that the software was used during solving the appropriate integer graph can not be colored with the given number of colors. programming. Throughout the testing the following, typ- Therefore the algorithm had to be able to bring in new col- ical basic regulations were used: ors. The DSATUR algorithm gave a proper solution for the coloring, whereabout the article of Brelaz [2] written in 1979 writes in details. In this case the algorithm was modified by • The minimum rest time between two shifts is 12 hours. using the set of employees estimated during the graph col- • The worktime should be less than 48 hours during one oring. Consequently a k-coloring was made in a way where week. adding a new color is possible. The graph coloring gives an initial solution which meet the given constraints and as- • Employees should have at least 4 free days in one signe every shift to a driver. As we know the duration of month. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 11 • The number of the maximum consecutive workdays is 6. CONCLUSION AND FUTURE WORKS 6. The important aspects of the crew rostering in the schedul- ing of bus drivers were introduced above. We proposed a two-phase graph coloring method for crew rostering. In the In our case every employee had a contract with 8 hours of first step, a graph was built and colored, while in the second worktime. Artificially generated and real life inputs were step, the graph was recolored with the tabu search method used in the test with one month planning period. The inputs by our algorithm. In the majority of test cases, the algo- being used during the test had a worktime between 7 and rithm has reached the theoretical lower bound. Our method 9 hours, and the duty time between 6 and 10 hours. While has been tested with artificially generated and real inputs. solving the problem, the objective function is divided into For moderate size problems, the results of the new algorithm two parts. The first is the set of people which is estimated have been compared to the solutions of the appropriate in- in the initial phase, while the minimization of the overtime teger programming model. The heuristic produced a satis- is considered during the tabu search. The tables below show factory running time, reached the lower bount in most cases the under- and overtime by optimal employee number and and returned a feasible solution in all cases. In the future the the maximum amount of daily shifts being 25. The results running time could be improved with paralell programming. by the input size being 50 and by real input are listed below In addition, the ability of handling drivers with different as well. kind of contracts with the heuristic would be desired. Table 1: Test results on generated inputs 7. ACKNOWLEDGEMENT Lower Time Cost Time Cost This work was partially supported by the European Union Input bound IP IP heur heur and the European Social Fund through project Supercom- gen25 1 4807 18.77 4900 3.2 4807 puter, the national virtual lab (grant no.: TAMOP-4.2.2.C- gen25 2 4658 495.2 4691 3.6 4658 11/1/KONV-2012-0010). gen25 3 1825 154.8 1855 4.1 1825 I would also like to thank the Tisza Volan for the test data gen25 4 4371 22.96 4457 3.1 4371 and the technical consultation. gen25 5 4815 18.11 4911 5.4 4815 8. REFERENCES [1] J. Bartholdi. A guaranteed-accuracy round-off The stopping criteria of the algorithm was given to stop algorithm for cyclic scheduling and set covering. when it could not find a better solution after 200 iteration, Operation Research, 29:501–510, 1981. or there were only drivers with under- or overtime left. [2] D. Brélaz. New methods to color vertices of a graph. 22(4):251–256, Apr. 1979. Table 2: Test results on generated inputs [3] S. R. D. Sitompul. Nurse scheduling models: a Lower Time Cost Time Cost state-of-the-art review. J. Soc. Health Syst., Input bound IP IP heur heur 7:441–499, 2004. gen50 1 3359 875.1 3465 5.4 3359 [4] G. Dantzig. A comment on edies traffic delays at toll gen50 2 2581 844.8 2581 15.2 2581 booths. Operations Research, 2:339–341, 1954. gen50 3 3000 1102.34 3047 10.4 3000 [5] K. M. M. H. Dowling, D. and D. Sier. Staff rostering gen50 4 1671 1252.78 1696 4.8 1671 at a large international airport. Annals of Operations gen50 5 1925 1529.56 1960 10.3 1925 Research, 72:125–147, 1997. [6] M. Elshafei and H. K. Alfares. A dynamic programming algorithm for days-off scheduling with The real examples belonged to the database of the Tisza sequence dependent labor costs. J. of Scheduling, Volan and consisted the local trips, shifts and the drivers in 11(2):85–93, Apr. 2008. Szeged. The real input consisted approximately 3000 shifts and the employees had a maximum number of 200 during [7] M. Gamache, A. Hertz, and J. O. Ouellet. A graph the solution. The results of the real input were tested with coloring model for a feasibility problem in monthly appropriate integer programming model, but there was no crew scheduling with preferential bidding. Comput. solution by the 2 % gap. The best results can be seen in the Oper. Res., 34(8):2384–2395, Aug. 2007. following table. [8] M. R. Garey and D. S. Johnson. Computers and Intractability: A Guide to the Theory of NP-Completeness. W. H. Freeman & Co., New York, Table 3: Test results on real life input NY, USA, 1979. Lower Time Cost Time Cost [9] F. Glover and M. Laguna. Tabu Search. Kluwer Input bound IP IP heur heur Academic Publishers, Norwell, MA, USA, 1997. volan 4897 21418.56 5012.87 9.87 4897 [10] D. Marx. Graph colouring problems and their applications in scheduling, 2003. The test results show that the task could be solved with [11] K. Nurmi, J. Kyngäs, and G. Post. Driver rostering the set of estimated number of employees in every case. We for bus transit companies. Engineering Letters, obtained that our algorithm is able to handle relatively large 19(2):125–132, December 2011. inputs, and in the majority of the test cases, it has reached [12] S. O. T. R. T. Grossman, D. Samuelson. Call centers. the theoretical lower bound with producing a satisfactory Technical Report, Haskayne Sch. Business,Univ. running time, and it produced a feasible solution in all cases. Calgary, 1999. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 12 A New Heuristic Approach for the Vehicle and Driver Scheduling Problems Viktor Árgilán Balázs Dávid University of Szeged University of Szeged Árpád tér 2. Boldogasszony sgt. 6. Szeged, Hungary Szeged, Hungary gilan@jgypk.u-szeged.hu davidb@jgypk.u-szeged.hu ABSTRACT out. Both problems are NP-hard, and literature usually of- Both vehicle and driver scheduling, the two main problems fers two different approaches for their solution: sequential in the optimization of public transport, are NP-hard. Be- or combined. The operational costs of the company come cause of this, finding an exact solution for a real-life instance from the combined costs of the above two problems, which is usually not possible. As a result, research of this field is leads to a complex optimization task. concentrated on developing efficient heuristic solution meth- ods for these problems. Balogh and Békési give a method to The aim in both cases is to create schedules where are trips create driver schedules based on an existing vehicle schedule. are executed exactly once. The sequential approach creates However, the running time of their algorithm is large. In this a feasible vehicle schedule as a first step, and then uses this paper, we examine the possibility of a heuristic acceleration to as a basis for driver schedules, if it is possible. The com-technique, which is tested both on randomly generated and bined method aims to solve both problems in one step. This real data. would mean that the combined method can provide an op- timal solution, but the size of the problem is usually too Categories and Subject Descriptors big to handle. On the other hand, while the solution given by the sequential method most likely will not be optimal, H.4 [Information Systems Applications]: Miscellaneous; the problem can be handled more easily in two steps. Both D.2.8 [Software Engineering]: Metrics—complexity mea- the vehicle and driver scheduling problems are NP-hard [2], sures, performance measures which makes it difficult to solve real problems with a large number of trips. This lead to a research into different heuris-General Terms tics for these problems, for example in [9, 11]. Using the idea Theory found in [7], Balogh and Békési propose a sequential solution algorithm for the problem in [1]. We introduce a heuristic Keywords acceleration technique that can speed up the running time of their algorithm. Vehicle schedule, driver schedule, sequential approach, com- bined approach In this paper, we first give a overview of the vehicle and driver scheduling problems, then briefly overview their se- 1. INTRODUCTION quential application. After this, we introduce our heuristic Cost reduction is important for public transport companies, method, which we apply to the solution methodology pro- and this can mostly be achieved by reducing operational posed in [1]. We present the efficiency of our heuristic both costs. These costs are derived from the long-term schedules on random and real instances. created by the company. These are then executed on a daily basis. The optimization of these schedules is a difficult prob-2. THE VEHICLE AND DRIVER SCHEDUL- lem, which is based on the pre-determined timetable of the ING PROBLEM company. The timetable consists of trips, which are carried out on a fixed route. There are two main problems connected Related works for the vehicle scheduling problem concern to the trips: vehicle and driver scheduling. Vehicle schedul- two different groups, depending on the types of available ve- ing gives the duties that different vehicles have to execute, hicles and the geographical location of the depots. The sin- while driver scheduling defines the shifts that driver carry gle depot vehicle scheduling problem (SDVSP) has a single vehicle type, and every vehicle starts at the same geographi- cal location. This problem can be solved in polynomial time. If there are more vehicle types or starting geographical lo- cations, then the problem has at least 2 depots. This case is called a multiple depot vehicle scheduling problem, which is NP-hard [2]. Vehicle scheduling problems arising in the real world usually belong to this group. The vehicle scheduling problem assigns vehicles to the timetabled trips, so that the following criteria are met: StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 13 • Every trips has to be executed exactly once. programming model for this problem, which is originally pre- sented in [1]: • Every trip is compatible with the depot and vehicle type of its assigned vehicle. cdkxdk (1) • Each trip ti and tj assigned to the same vehicle have d∈D k∈Kd to be compatible. Two trips ti and tj are compatible, if they can be served by the same vehicle type, and xd there is enough time between the arrival time of t k ≥ 1f or∀t ∈ T (2) i and d∈D the departure time of t k∈Kd(t) j for the vehicle to travel from the arrival location of ti to the departure location of tj . Such a vehicle journey (without executing a timetabled xdk ∈ 0, 1, d ∈ D, k ∈ Kd (3) trip) is called a deadhead trip. • The problem wants to minimize the arising cost, which where is usually the combination of two main components: – a daily cost for each vehicle, if it leaves its depot • T is the set of timetabled trips, to execute a trip, • D is the set of depots, – and a travel cost for each vehicle, which is pro- portional to the travelled distance. • Kd is the set of possible driver schedules from depot d, Many different mathematical models can be given for the • Kd(t) is the set of driver schedules covering trip t from vehicle scheduling problem. For an overview of the possible depot d, representations, refer to [3]. • cdk is the cost of schedule k from depot d. Driver scheduling is also a main problem in public transit, because personell give a high percent of the operational costs of a company. The main problem of is model is that the number of pos- sible driver schedules is high, and generating all of them is The problem considers a set of driver tasks (which are mainly usually not possible. To address this problem, the column the timetabled trips), and the aim is to return driver shifts, generation method in [4] is applied, generating new columns which represent a sequence of the tasks. Similarly to vehicle using the time space network model in [8]. However, the scheduling, every task has be assigned to a shift exactly once, running time of this method is still large. and the tasks belonging to the same shift have to be com- patible. The most important constraints for driver schedules 4. THE PROPOSED HEURISTIC are defined by the European Union, but companies usually As we mentioned at the end of the last section, running define their own set of rules also. The most important of time of the above method is still an issue: it can be as big as these constraints are: several days. This running time might not be acceptable for a solution method applied in practice. Our aim is to develop • maximum working hours, a method that decreases the size of the problem itself, and results in a model with a smaller number of variables, which • total length of the daily break(s), is easier to solve. • and the number and length of the breaks. Mathematical formulation of the problem is usually given by set partitioning or a set covering relaxation. These problems are known to be NP-hard [6]. For an overview of the existing models and methods, refer to [5]. 3. SOLUTION WITH COLUMN GENERA- TION Balogh and Békési give a sequential solution method for ve- hicle and driver scheduling in [1]. Their methodology is strongly based on the ideas of Gintner et al. [7] and Steinzen et al. [10]. This method considers an arbitrarily given vehi-Figure 1: Traffic load of a sample problem cle schedule, and uses this as an initial step to construct the driver schedule. The decrease in size is achieved by the following method: consider the timeline of the daily traffic load as our plan- As we mentioned earlier, set covering is one of the usual ning period, which is defined by the timetabled trips. This approaches to model this problem. We give a 0-1 integer planning period starts with the earliest trip departure time, StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 14 As one of the most important driver rules is the assignment of breaks, we also have to check if there is enough idle time between the tasks of a vehicle duty to assign a break. If there is such a possibility, then dbt = 1 for a merged trip t. These newly constructed trips will be added to a set M . The heuristic continues sequentially with the sets Si (1 < i ≤ k). For each step i, the vehicle scheduling is carried out on the set Si ∪M . After the scheduling is completed, a newly merged trip t defined by a driver duty will have dbt = dbt, t where t are the trips of the duty. If there is enough time for and additional driver break in the new duty, dbt is further Figure 2: An interval division of the problem increased by 1. Set M will be then updated to only contain the trips defined by this new duty. and ends with the latest trip arrival time. An example for We explain the algorithm on the example given in Figure 3, this can be seen in Figure 1. where the set of trips is given for the problem over a timeline. We divide the planning period into smaller, fixed-length in- tervals. The heuristic will consider these sequentially, starting with the earliest interval. An example for such a division can be seen in Figure 2. The aim of the heuristic is to decrease the number of trips in the input, which is achieved by merging more trips together, which are then considered as a single trip. First, as we mentioned earlier, we divide our planning period into smaller intervals of equal size t. The number k of dividing points Figure 3: Input of the problem needed for this can be calculated by Suppose that the first dividing point for the heuristic will be at k1. We will only consider trips starting before k1, as (max{arrival time k = trip }−min{departure timetrip }) can be seen in Figure 4. t Using these dividing points, we partition the trips of the in- put into sets Si (1 ≤ i ≤ k). A set Si contains the trips that depart in the time interval defined by its dividing points. The set Sk will contain trips starting after the last dividing point k. It is possible, that there are only arriving trips, and no trips are departing after time k, which will result in Sk = ∅. Naturally, Si Sj = ∅ if i = j Figure 4: First iteration of the algorithm and k The vehicle scheduling problem is solved for these trips, and Si = S some of them are merged together, forming the new trips of i=1 set M . The merged trips that come as a result of the first Our heuristic aims to solve the vehicle scheduling problem phase can be seen in Figure 5 sequentially for every interval Si. An arbitrary solution method can be used for this. As the newly merged trips act more like the pieces of a shift, each trip t is assigned a number dbt, which represents the number of driver breaks that can be carried out between its sub-trips. Naturally, this value will be 0 for normal trips, but it can be positive for merged trips. The vehicle scheduling problem is solved for the trips belong- ing to set S1. This will result in a series of vehicle duties, each such duty consisting of at least one trip. The trips of a duty will be merged together to form a new trip. The Figure 5: Merged trips after the first iteration departure time and location of this new trip will correspond to that of the first trip of the duty, and the arrival time and Using these newly created trips, we now advance to the sec-location will be defined by the last trip of the duty. ond iteration, which has the dividing point k2. This iteration StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 15 considers the original trips departing between k Algorithm 1 Heuristic Size Reduction. 1 and k2, to- gether with the newly merged trips of set M . This can be 1: procedure SizeReduce(S, t, B, p) seen in Figure 6. 2: M=∅(max{arrival time 3: k = trip }−min{departure timetrip }) t 4: for i = 0 to k do 5: if i ∗ t < arrival timetrip ≤ (i + 1 ∗ t) then 6: dbtrip = 0 7: trip ⇒ Si 8: end if 9: end for 10: for i = 0 to k do 11: X = M Si 12: M = solve vehicle scheduling for X Figure 6: Second iteration of the algorithm 13: for j = 1 to m’ do 14: Merge tasks on duty j to a single trip 15: Calculate db The vehicle scheduling problem is once again solved for the j 16: end for trips, and the resulting vehicle duties are merged into single 17: M = M trips. The results of this iteration can be seen in Figure 7. 18: end for 19: end procedure Table 1 presents our results, while we present their solutions in Table 2 as a comparison. The driver rules we considered were the ones in practice at the transportation company. Table 2: Test results from [1] Problem Trips Vehicles Drivers Running time Figure 7: Merged trips after the second iteration #1 100 5 9 81min #2 900 38 62 362min The results in Figure 7 also present the consideration of #3 2711 103 224 6616min driver breaks. One of the merged trips has enough time for a break (colored green in the figure). As a result, the It can be seen in the last column that the running time of possibility of the break is noted by the assigned variable. the new method is smaller, even for large examples. Our new heuristic is faster for every instance than the old one. The pseudo-code of the algorithm can be seen in Algorithm However, it can also be seen in the third and fourth column 1. The algorithm first calculates the dividing points and as- that the number of vehicles and drivers given by our method signs trips to every set Si, then sequentially solves the vehi-are higher in every instance, which results in a larger cost. cle scheduling problem for these sets. The possible number As we expected, the heuristic managed to speed up the so- of driver breaks is calculated using the method described lution process significantly, but it produces solutions with a above. higher cost. 5. COMPUTATIONAL RESULTS 6. CONCLUSION AND FUTURE WORKS We tested the above size reduction heuristic on both ran- This is a heuristic acceleration of an earlier method. (of domly generated and real-life data provided by Tisza Volán Balogh and Bekesi). It produces less running time and a bit Zrt, the bus company of the city Szeged, Hungary. The higher cost, but the gap between the value of our solution data we used was the same that Balogh and Békési used for and the optimal value is small. As a future work, we plan testing in [1]. The tests run on the following computer: to apply a new matching based method instead of greedy heuristic. And we plan to examine other heuristics as well. • Processor: Intel Core I7 X980 3.33 Ghz, 7. ACKNOWLEDGEMENT • Memory: 12Gb, This work was partially supported by the European Union • and the European Social Fund through project Supercom- Operation System: Microsoft Windows 7 Ultimate 64bit puter, the national virtual lab (grant no.: TAMOP-4.2.2.C- 11/1/KONV-2012-0010). Table 1: Test results of heuristic method 8. REFERENCES Problem Trips Vehicles Drivers Running time [1] J. Balogh and J. Békési. Driver scheduling for vehicle #1 100 5 9 69min schedules using a set covering approach: a case study. #2 900 39 65 264min Abstracts of the ICAI 2014 - 9th International #3 2711 105 231 4989min Conference on Applied Informatics, to appear. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 16 [2] A. A. Bertossi, P. Carraresi, and G. Gallo. On some matching problems arising in vehicle scheduling models. Networks, 17(3):271–281, 1987. [3] S. Bunte and N. Kliewer. An overview on vehicle scheduling models. Journal of Public Transport, 1(4):299–317, 2009. [4] M. Desrochers and F. Soumis. A column generation approach to the urban transit crew scheduling problem. Transportation Science, 23(1):1–13, 1989. [5] A. T. Ernst, H. Jiang, M. Krishnamoorthy, and D. Sier. Staff scheduling and rostering: A review of applications, methods and models. European journal of operational research, 153(1):3–27, 2004. [6] M. Fischetti, A. Lodi, S. Martello, and P. Toth. A polyhedral approach to simplified crew scheduling and vehicle scheduling problems. Management Science, 47(6):833–850, 2001. [7] V. Gintner, N. Kliewer, and L. Suhl. A crew scheduling approach for public transit enhanced with aspects from vehicle scheduling. In Computer-aided Systems in Public Transport, pages 25–42. Springer, 2008. [8] N. Kliewer, T. Mellouli, and L. Suhl. A time–space network based exact optimization model for multi-depot bus scheduling. European journal of operational research, 175(3):1616–1627, 2006. [9] A.-S. Pepin, G. Desaulniers, A. Hertz, and D. Huisman. A comparison of five heuristics for the multiple depot vehicle scheduling problem. Journal of Scheduling, 12(1):17–30, 2009. [10] I. Steinzen, V. Gintner, L. Suhl, and N. Kliewer. A time-space network approach for the integrated vehicle-and crew-scheduling problem with multiple depots. Transportation Science, 44(3):367–382, 2010. [11] A. Tóth and M. Krész. An efficient solution approach for real-world driver scheduling problems in urban bus transportation. Central European Journal of Operations Research, 21(1):75–94, 2013. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 17 StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 18 Multi-population Firefly Algorithm Jani Dugonik Iztok Fister Faculty of Electrical Engineering and Computer Faculty of Electrical Engineering and Computer Science Science University of Maribor University of Maribor 2000 Maribor, Slovenia 2000 Maribor, Slovenia jani.dugonik@um.si iztok.fister@um.si ABSTRACT thus learn the success of problem-solving from nature and This paper proposes a meta-heuristic Multi-Population Fire- develop nature-inspired heuristic and/or meta-heuristic al- fly Algorithm (MPFA) for single-modal optimization using gorithms in order to solve optimization problems with which two multi-population models, i.e., one is based on the island developers are confronted today. Two sources from the na- model while the other on the mainland-island model. The ture have particularly inspired developers of new optimiza- unique characteristics of each sub-population is evolved in- tion algorithms, i.e., Darwinian evolution and the behavior dependently and the diversity of the entire population is ef- of social living insects (e.g., ants, bees, termites, etc.) and fectively increased. Sub-populations communicate with each other creatures (e.g., birds, dolphins, fireflies, etc.). As a other to exchange information in order to expand the search result, two main classes of nature-inspired algorithms exist range of the entire population. In line with this, each sub- nowadays, i.e., evolutionary algorithms (EA) [3] and swarm population explores a specific part of the search space and intelligence (SI)-based algorithms [1]. While the former al- contributes its part for exploring the global search space. ready came in mature years, the latter has experienced rapid The main goal of this paper was to analyze the performance development. Almost every day, we are witnessing the birth between MPFA and the original Firefly Algorithm (FA). Ex- of a new SI-based algorithm. periments were performed on a CEC 2014 benchmark suite consisting of 16 single-objective functions and the obtained One of the younger members of the SI-based algorithms is results show improvements in most of them. the Firefly Algorithm (FA) as proposed by Yang in [11]. FA is inspired by a chemical phenomenon bioluminiscence Keywords needed by natural fireflies to find their prey, on the one hand, and to attract their mating partners, on the other swarm intelligence, island model, multi-population, firefly hand. This algorithm belongs to a class of population-based algorithm algorithms [5, 4, 11, 12]. Population-based approaches main- tain and improve multiple candidate solutions, often using 1. INTRODUCTION population characteristics to guide the search. Two major An optimization problem is defined as a quadruple OP = components of any population-based search algorithms are I, S, f, goal , where I denotes a set of all input instances exploitation and exploration. Exploitation refers to search- x ∈ I in the form of x = {x ing within neighborhoods for the best solutions and ensures 1, x2, ..., xD } where D is the dimensionality of the problem, S(x) a set of all feasible so- that the solutions can converge into optimality, while the lutions y = S(x), f is the objective function estimating exploration uses randomization in order to avoid the solu- the feasible solution, and goal determines an optimal cri- tions being trapped within a local optima and while at the teria that can be either the minimum or maximum value same time increasing the diversity of the solutions. A good of the objective function. A task of the optimization algo- combination of these two components may usually ensure rithm is to find the value of y∗ that minimizes or maximizes that the global optimality is achieved [11]. (depending on the goal ) the value of the objective function f (y). The domain values of input variables x One of the possible ways of how to improve the exploration i ∈ [lb i, ub i] are limited by their lower lb and exploitation in the original FA algorithm can be split- i and upper ubi bounds. ting the FA population into more sub-populations (so-called Nature has evolved over millions of years and has found per- multi-populations). For instance, the authors in [10] pre- fect solutions to almost all encountered problems. We can sented a multi-population FA for correlated data routing within underwater wireless sensor networks. They designed three kinds of fireflies and their coordination rules in or- der to improve the adaptabilities of building, selecting, and optimizing a routing path. Groups are represented as sub- populations, where each sub-population conducts its own optimization in order to improve the convergence speed and solution precision of the algorithm. The author in [13] ana- lyzed the ability of a multi-population differential evolution to locate all optima of a multi-modal function. The explo-StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 19 ration was ensured by the controled initialization of sub- Algorithm 1 Firefly Algorithm populations while a particular differential evolution algo- 1: Objective function f (x), x = (x1, ..., xD)T rithm ensured the exploitation. Sub-populations were com- 2: Generate initial population of fireflies xi (i=1,2,...,N p) municating via archive where all located optima were stored. 3: Light intensity Ii at xi is determined by f (xi) The authors in [8] used an evolutionary algorithm on punc- 4: Define light absorption coefficient γ tuated equilibria. The theory of punctuated equilibria calls 5: while (t < Gmax) do for the population to be split into several sub-populations. 6: for i=1 to n f iref lies do These sub-populations have isolated evolution (computa- 7: for j=1 to n f iref lies do tion) and scattered with migration (communication). 8: if (Ij > Ii) then 9: Move firefly i towards firefly j Eq. (3) This paper aimed to evaluate whether it is possible to out- 10: end if perform the performance of the original FA algorithm by 11: Evaluate new solution and update light in- splitting its population into more sub-populations. The pro- tensity Eq. (1) posed multi-population FA (MPFA) supports two multi- 12: end for population models, i.e., Island [13] and Mainland-Island [9]. 13: Rank fireflies and find the current best In these multi-population models, sub-populations evolve in- 14: end for dependently, thus the unique characteristics of each sub- 15: Post-process results and visualization population can be effectively maintained, and the diver- 16: end while sity of the entire population is effectively increased. Sub- populations communicate with each other by exchanging in- formation in order to expand the search range of the entire in Algorithm 1. Light intensity of a firefly is defined as: population. The search technique based on a population has proved to have good ability regarding global searching and I(r) = I0 · e−γ·r2 (1) can find a set of solutions in one-shot operation. The pro- posed multi-population FAs were compared with the original where I FA on single-objective CEC 2014 benchmark functions [7]. 0 is the original light intensity at the location of r = 0, γ is the light absorption coefficient and r is the distance between two fireflies. The distance between any two fireflies The remainder of this paper is organized as follows. In Sec- i and j at x tion 2 the original FA will be presented. In Section 3 a i and xj can be expressed as Cartesian distance r multi-population FA with two multi-population models are ij = ||xi − xj ||. As firefly attractiveness is proportional to the light intensity, we can define the attractiveness of a presented in detail. Section 4 presents experiments, where firefly using the following equation: the number of tests were performed in order to compare the proposed approach with the original FA. The paper is con- β(r) = β0e−γ·r2 (2) cluded with Section 5, where our opinion on the obtained results is given. where β0 is their attractiveness at r = 0. Firefly i that is 2. THE ORIGINAL FIREFLY ALGORITHM attracted to another more attractive firefly j is determined by: The Firefly Algorithm (FA) [11] has two fundamental fac- tors: light intensity and attractiveness. Light intensity I xi = β0e−γ·r2 · (xj − xi) + α · εi (3) reflects the firefly location and determines its direction of movement, while the degree of attractiveness determines the which is randomized with the vector of random variable ε distance that a firefly has moved. Both factors are con- i, being drawn from a Gaussian distribution, and step factor stantly updated in order to achieve the objective of the op- α ∈ [0, 1]. timization. For simplicity, the author in [11] used the following three 3. THE MULTI-POPULATION FIREFLY AL- idealized rules: GORITHM The multi-population firefly algorithm (MPFA) can be sum- marized in the pseudo-code as shown in Algorithm 2. MPFA • All fireflies are unisex so that one firefly will be at- will consider there to be an overall population P of Np tracted to other fireflies regardless of their sex. fireflies (individuals) that is split into N sub-populations • Attractiveness is proportional to their brightness and P1, P2, ...PN . Each sub-population has Nsp individuals and for any two flashing fireflies, the dimmer one will move the number of N sub-populations is calculated with the fol- towards the brighter one. They both decrease as their lowing equation: distance increases. If there is no brighter one, it will Np N = (4) move randomly. Nsp • The brightness of a firefly is affected or determined by the landscape of the objective function. For sub-populations to communicate with each other, the magnitude and frequency of that communication are neces- sary. These two parameters determine the amount of iso- Based on these three rules, the basic steps of the firefly al- lation and interaction between sub-populations. Periods of gorithm (FA) can be summarized as the pseudo-code shown isolated evolution are referred to as epoch, with migration StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 20 Algorithm 2 Multi-Population Firefly Algorithm 1: Objective function f (x), x = (x1, ..., xD)T 2: Calculate number of sub-populations (N ) 3: for all n ∈ N do 4: Generate initial sub-population Pn of fireflies xi (i=1,2,...,Np) 5: end for 6: Light intensity Ii at xi is determined by f (xi) 7: Define light absorption coefficient γ 8: for e=1 to Gmax do epoch Figure 1: Island Model 9: for all n ∈ N do 10: for g=1 to epoch do 11: for i=1 to n f iref lies do Algorithm 3 Multi-Population Firefly Algorithm with Is- 12: for j=1 to n f iref lies do land Model - Migration 13: if (Ij > Ii) then 14: Move firefly i towards firefly j 1: Get the number of migrants Nm to migrate 15: end if 2: for i = 1 to N do 16: Evaluate new solutions 3: Choose Nm individuals from Pi that are mutually 17: Update light intensity different and save them to the matrix M 18: end for 4: end for 19: Rank fireflies and find the current best 5: Mutually exchange individuals that are defined in ma- 20: end for trix M 21: end for 22: end for 23: Migrate fireflies of mainland and islands, where mainland and islands are re- 24: end for ferred to as sub-populations. When each sub-population has 25: Find the best firefly from all sub-populations executed a sequential FA for Gi generations, individuals are 26: Post-process results and visualization migrated from sub-populations P2, ..., PN to sub-population P1, as shown in Algorithm 4. Let us assume two sub- populations P1 and P2 with Nsp = 10 and Nm = 20% are de- occurring at the end of each epoch except the last. The fined. Then, two individuals per sub-population P 2, ..., PN length of the epoch determines the frequency of interac- are moved to sub-population P1, as shown in Figure 2. At tion and is usually specified by a number of generations the end of migration, sub-population P1 is sorted according (epoch) that P to the fitness values of migrated individuals. In order to keep n evolves in isolation. During the epoch, each sub-population executes a sequential FA for epoch indepen- the size of sub-population P1 the same, the top Nsp individ- dently. At the end of each epoch, individuals are migrated uals are retained, while the others are discarded. After the between sub-populations. There are many various migra- algorithm has reached the terminating criteria, the best in- tion strategies in multi-population models. The following dividual was taken from the sub-population P1 (mainland). two models are described and used in this paper: island and mainland-island model. Algorithm 4 Multi-Population Firefly Algorithm with Mainland-Island Model - Migration 3.1 Island Model 1: Get the number of migrants Nm to migrate 2: for i = 2 to N do The island Model Firefly Algorithm (MPFA-In, where n de- 3: Choose N termines the number of sub-populations) consist of islands, m individuals from Pi that are mutually different and save them to the matrix M where islands are referred to as sub-populations. When each 4: end for sub-population is executed a sequential FA for epoch gener- 5: Copy chosen individuals from sub-populations ations, individuals are migrated between sub-populations, P as shown in Algorithm 3. The magnitude of the commu- 2, ..., PN to sub-population P1 6: Sort individuals in P nication is defined, for instance, as N 1 by light intensity m = 25%. Then, 7: Keep top Nsp individuals, remove others, so the size of Nm percent of migrants are chosen for each sub-population the sub-population P which were exchanged with other sub-populations, as shown 1 remains the same in Figure 1. Let us assume two sub-populations P1 and P2 with Nsp = 10 and Nm = 20% are defined. Then, two 4. EXPERIMENTS AND RESULTS individuals from sub-population P1 are exchanged with two The goal of this experimental work was to show that MPFA individuals in sub-population P2. Thus, the sizes of the sub- can outperform the results of the original FA algorithm. In populations remain the same. After the algorithm reaches our experiments, the results of the original FA were com- the termination criteria, the best individual is taken from pared with the results of the following MPFA: MPFA-I2 all sub-populations. and MPFA-I4 (i.e., MPFA with island model using two or four sub-populations), and MPFA-M2 and MPFA-M4 (i.e., 3.2 Mainland-Island Model MPFA with mainland-island model using two or four sub- The Mainland-Island Model Firefly Algorithm (MPFA-Mn, populations). Additionally, the following three population where n determines the number of sub-populations) consist models were used during tests, i.e., small with an original StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 21 Table 1: Comparison between FA algorithms for population model Np=100 and D=10 Func. FA MPFA-I2 MPFA-I4 MPFA-M2 MPFA-M4 1 1.0243e+06 ± 2.0147e+06 7.6582e+05 ± 3.4941e+06 6.4129+e05 ± 2.3202e+06 7.6582e+05 ± 3.4929e+06 5.0199e+05 ± 9.3012e+05 2 1.3862e+04 ± 8.2668e+03 6.1151e+03 ± 4.8489e+03 5.4070+e03 ± 4.0991e+03 6.1151e+03 ± 4.8489e+03 6.5121e+03 ± 8.0172e+03 3 2.0877e+04 ± 1.9394e+04 2.2167e+04 ± 2.0586e+04 1.6543e+04 ± 1.3981e+04 2.2167e+04 ± 2.0586e+04 2.3253e+04 ± 2.6483e+04 4 7.7409e+00 ± 1.6024e+01 6.9890e+00 ± 1.1157e+01 7.0024e+00 ± 2.9768e+01 7.1320e+00 ± 1.0983e+01 6.8717e+00 ± 8.8252e+00 5 2.0107e+01 ± 0.0704e+00 2.0059e+01 ± 0.0405e+00 2.0035e+01 ± 0.0260e+00 2.0059e+01 ± 0.0405e+00 2.0044e+01 ± 0.0380e+00 6 8.4299e+00 ± 4.0688e+00 8.9432e+00 ± 2.5879e+00 8.3829e+00 ± 2.6099e+00 8.9432e+00 ± 2.5879e+00 9.6953e+00 ± 3.2605e+00 7 5.4650e+00 ± 6.2455e+00 9.9181e+00 ± 1.0126e+01 5.8841e+00 ± 4.6803e+00 9.9181e+00 ± 1.0126e+01 5.3097e+00 ± 5.9105e+00 8 1.6925e+01 ± 1.5711e+01 2.3883e+01 ± 2.3057e+01 2.1892e+01 ± 1.4542e+01 2.3883e+01 ± 2.3057e+01 3.0847e+01 ± 2.5583e+01 9 1.6924e+01 ± 1.2980e+01 1.8907e+01 ± 2.0379e+01 1.9903e+01 ± 1.6285e+01 1.8907e+01 ± 2.0379e+01 3.0847e+01 ± 2.7831e+01 10 1.1733e+03 ± 9.0738e+02 1.3110e+03 ± 1.1059e+03 1.0351e+03 ± 8.9937e+02 1.3110e+03 ± 1.1157e+03 1.2901e+03 ± 1.4476e+03 11 1.1434e+03 ± 9.9968e+02 1.2193e+03 ± 1.0618e+03 9.1012e+02 ± 8.3705e+02 1.2193e+03 ± 1.0618e+03 1.3165e+03 ± 1.1339e+03 12 0.4707e+00 ± 0.9716e+00 0.3272e+00 ± 1.0509e+00 0.1942e+00 ± 0.4033e+00 0.3272e+00 ± 1.0509e+00 0.4343e+00 ± 1.7646e+00 13 0.3973e+00 ± 0.3133e+00 0.3949e+00 ± 0.3414e+00 0.3137e+00 ± 0.2133e+00 0.3949e+00 ± 0.3414e+00 0.3578e+00 ± 0.3217e+00 14 0.3618e+00 ± 0.2481e+00 0.3578e+00 ± 0.2386e+00 0.3260e+00 ± 0.1413e+00 0.3578e+00 ± 0.2386e+00 0.3378e+00 ± 0.2627e+00 15 1.3559e+01 ± 1.5484e+01 1.7627e+01 ± 1.5247e+01 2.1337e+01 ± 1.8450e+01 1.7627e+01 ± 1.5247e+01 2.8024e+01 ± 2.6350e+01 16 3.7235e+00 ± 0.5893e+00 3.8952e+00 ± 0.7525e+00 3.7126e+00 ± 0.7731e+00 3.8952e+00 ± 0.7525e+00 4.0390e+00 ± 0.7963e+00 tion. Note that error values smaller than 10−8 were taken as zero. In order to limit a search space, each problem variable can capture the value from the range xi ∈ [−100, 100]D, where values -100 and 100 represent its upper and lower bounds. The dimensionality of all problems was limited to D = 10. The results of the mentioned FA algorithms are illustrated in Table 1. The results for all population models (i.e., small, medium, large) were obtained. Small population model (Np = Figure 2: Mainland-Island Model 100) gave the best results, and due to the paper’s length lim- itation only these results are presented in this table. In line with this, the original FA algorithm is compared with the MPFA-I2 and MPFA-M2 using two sub-populations of size population size of Np = 100, medium with Np = 200 and Nsp = 50, and MPFA-I4 and MPFA-M4 using four sub- large with Np = 400. The original population size was di- populations of size Nsp = 25. The table presents mean and vided between sub-population according to Eq. (4). The standard deviation over 51 independent runs for each algo- same number of generations Gmax = 1, 000 was used for rithm. The results in Table 1 show that MPFA-I2 as well each sub-population for each algorithm. as MPFA-M2 outperformed the original FA on 7 out of 16 test functions, MPFA-M4 on 8 out of 16 test functions and The FA algorithms used the following parameter settings. MPFA-I4 on 12 out of 16 test functions. The best results The maximum number of evaluations was set as MAX FEs = were obtained with MPFA-I4 which outperformed the other 10, 000 · D. The randomized factor was fixed at α = 0.5, the MPFAs and the original FA on 10 out of 16 functions. lights absorption at γ = 1, and the attractiveness at the be- ginning β0 = 1. In each generation the randomized factor α In order to evaluate the quality of the results statistically, was updated by the following equation: α = α·(1−δ), where 1 Friedman tests [6] were conducted that compare the average δ = 1.0 − ( 104 ) G [5]. For the MPFA algorithms, some ad- 0.9 ranks of the compared algorithms. Thus, a null-hypothesis is ditional parameters were used, like the number of epochs as placed that states: two algorithms are equivalent and there- epoch = 100, and migration probability Nm = 25%. Tests fore, their ranks should be equal. When the null-hypothesis were conducted on all three population models using five is rejected, the Bonferroni-Dunn test [2] is performed. In FA algorithms, i.e., FA, MPFA-I2, MPFA-I4, MPFA-M2, this test, the critical difference is calculated between the and MPFA-M4. In summary, 15 tests were performed, in average ranks of those two algorithms. If the statistical dif- which 51 independent runs were performed. ference is higher than the critical difference, the algorithms are significantly different. All algorithms were tested on the 16 single-objective uni- modal and simple multi-modal CEC 2014 benchmark func- Three Friedman tests were performed regarding data ob- tions [7]. For uni-modal functions the convexity guarantees tained by optimizing 16 functions of three different popula- that the final optimal solution is also the global optimum. tion sizes according to five measures. As a result, each algo- The global maximum was measured according to an error rithm during the tests (also the classifier) was compared with value ER. The error value for each function is calculated regard to the 48 functions x 5 measurements this means, by subtracting the value of global optima from the obtained 240 different variables. The tests were conducted at a sig- value according to the following equation: nificance level of 0.05. The results of the Friedman non- ER parametric test can be seen in Figure 3 that is divided into i = fi(x) − fi(x∗), [7] (5) three diagrams. Each diagram shows the ranks and confi- where i is the function number, fi(x) is the obtained value, dence intervals (critical differences) for the algorithms under and fi(x∗) = 100·i is the value of global optima for i-th func-StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 22 MPFA-M4 MPFA-M2 MPFA-I4 MPFA-I2 FA 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 (a) Small (Np = 100) (b) M edium (Np = 200) (c) Large (Np = 400) Figure 3: Results of the Friedman non-parametric test consideration with regard to the dimensions of the functions. [5] I. Fister, X.-S. Yang, I. Fister, and J. Brest. Memetic Note that the significant difference between the two algo- firefly algorithm for combinatorial optimization. In rithms is observed if their confidence intervals denoted as Proceedings of the Fifth International Conference on thickened lines in Fig. 3 do not overlap. Bioinspired Optimization Methods and their Applications - BIOMA 2012, pages 75–86, 2012. As can be seen from Fig. 3, the MPFA-I4 outperformed the [6] Milton Friedman. A comparison of alternative tests of results of the original FA as well as the other algorithms significance for the problem of m rankings. The Annals using all three observed population size model significantly. of Mathematical Statistics, 11:86–92, March 1940. The MPFA-M4 achieved the results that are significantly [7] J. J. Liang, B.-Y. Qu, and P. N. Suganthan. Problem worse than the results of the other FA algorithms. The Definitions and Evaluation Criteria for the CEC 2014 performances of the other three algorithms were comparable Special Session and Competition on Single Objective with each other. Real-Parameter Numerical Optimization. Technical report, Computational Intelligence Laboratory, 5. CONCLUSION Zhengzhou University, Zhengzhou China and In this paper, we proposed MPFA using two multi-population Technical Report, Nanyang Technological University, models, i.e., Island and Mainland-Island Models. The pro- Singapore, 2013. posed MPFAs were compared with the original FA algorithm [8] W. N. Martin, J. Lienig, and J. P. Cohoon. Island using three different population size models (i.e., small, medium, (migration) models: evolutionary algorithms based on large) by solving the CEC-14 benchmark function suite. Based punctuated equilibria. Handbook of Evolutionary on the obtained results, we can see that the most promising Computation, Oxford University Press, 1997. results were obtained by the MPFA-I4. This fact encourage [9] M. Slatkin and L. Voelm. FST in a Hierarchical Island us to continue with the experiments of multi-population FA Model. Genetics, 1991. in the future. [10] M. Xu and G. Liu. A Multipopulation Firefly Algorithm for Correlated Data Routing in Underwater The future work could be especially focused on the migra- Wireless Sensor Networks. International Journal of tion probability and dimension of the problem. Current mi- Distributed Sensor Networks, 2013. gration probability was fixed for all multi-population mod- [11] X.-S. Yang. Nature-Inspired Metaheuristic Algorithms. els, but the migration probability can be modified or even Luniver Press, 2008. adapted during the algorithm run. On the other hand, all [12] X.-S. Yang. Multiobjective firefly algorithm for the performed tests were done on small dimensions of the continuous optimization. Engineering with Computers, problem. Thus, the algorithm with few number of evalua- pages 175–184, 2013. tions and larger population sizes did not reach the migration [13] D. Zaharie. A Multipopulation Differential Evolution phase at all. With larger dimensions, the number of evalua- Algorithm for Multimodal Optimization. In tions would be increased and the multi-population strategies Proceedings of Mendel 2004 - 10th International could perform even better. Conference on Soft Computing, pages 17–22, 2004. 6. REFERENCES [1] C. Blum and D. Merkle. Swarm Intelligence. Springer-Verlag, Berlin, 2008. [2] Janez Demšar. Statistical comparisons of classifiers over multiple data sets. Journal of Machine Learning Research, 7:1–30, 2006. [3] A. Eiben and J. Smith. Introduction to Evolutionary Computing. Springer-Verlag, Berlin, 2003. [4] I. Fister, I. Fister Jr, X.-S. Yang, and J. Brest. A comprehensive review of firefly algorithms. Swarm and Evolutionary Computation, pages 34–46, 2013. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 23 StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 24 Monte Carlo Path Tracing with OpenCL Andrej Bukošek FRI, University of Ljubljana Večna pot 113, 1000 Ljubljana, Slovenia andrej.bukosek@gmail.com ABSTRACT 1.1 BRDF We introduce an interactive Monte Carlo path tracer that uses the OpenCL framework. A path tracer draws a photo- realistic image of a 3D scene by simulating physical effects of light. Interactivity enables the user to move around the n scene in real time, while OpenCL makes it possible to run ωi ωo the same code on either CPU or GPU architectures. The results presented here are a starting point for our further work on path tracing on heterogeneous architectures. Categories and Subject Descriptors I.3.7 [Computer Graphics]: Raytracing, Animation; D.1.3 [Software]: Parallel programming General Terms Figure 1: BRDF vectors. Algorithms, Performance Keywords The BRDF has two parameters: direction towards the light (ωi) and direction towards the camera (ωo), as shown in path tracing, OpenCL, rendering, computer graphics Figure 1. Based on these two directions and the internal parameters of the BRDF model, we can calculate the re-1. INTRODUCTION flectance of the object’s surface in the point hit by the ray of light. The path tracing algorithm draws a photorealistic image of a 3D scene by simulating the bouncing of light around the To obtain physically-correct results, the BRDF must satisfy scene in a physically-correct way [2]. Rendering consists of two properties: reciprocity and conservation of energy. A shooting rays of light from the camera into the scene, check- BRDF is reciprocal if it returns the same result if the light ing whether they intersect any object in their path, and re- and camera direction vectors are swapped. A BRDF con- cursively bouncing them off the intersected object based on serves energy if the amount of light reflected off an object physical laws. How the light is reflected off objects or passed is less than or equal to the amount of incident light. This through them depends on their material. We can describe means that a BRDF cannot generate light. this process using functions such as: BRDF (bidirectional re- flectance distribution function), BTDF (bidirectional trans- mittance distribution function), and BSSRDF (bidirectional 1.2 Lambertian BRDF surface scattering reflectance distribution function). The simplest BRDF is the Lambertian model, which doesn’t exist in nature (the closest materials available are unpro- BRDF describes how light reflects off opaque objects [2, 5] cessed wood and printer paper, but only if we look at it and defines the look of the material (e.g. whether it looks straight on and not from the sides). more metallic or more plastic). BTDF describes how light travels through a translucent object (e.g. glass), while BSS- RDF describes the scattering of light under the surface of ρ an object before it leaves the object (an example of such a BRDF (ωi, ωo) = (1) π material is human skin). Equation 1 shows the BRDF for this model. ρ represents the reflectance of the point hit by light. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 25 1.4 Numerical solution The path tracing algorithm shoots rays from the camera into the scene, bounces them around the objects within, and calculates the luminance of points hit along the way. A random sample of the hemisphere is needed when calcu- lating the contributions, since it’s impossible to check every conceivable direction. However, such an approach provably converges towards the right solution over time [3]. To improve the quality of generated images, supersampling is usu- ally used. The more samples we use, the less noise is present in the output image, however, the render time increases. 1.5 OpenCL Figure 2: Example of the Lambertian BRDF. OpenCL is a framework for executing programs written in a Figure 2 shows a sphere and walls made out of a Lambertian special C99-based language across heterogeneous platforms material. [4]. 1.3 Path tracing A platform can have many compute devices attached to a host processor, each of them can execute a kernel (a program We can calculate the luminance of each point in a 3D space written in OpenCL’s C-like language). Compute devices can using the global illumination equation (2). be standard processors (CPUs), graphics processing units (GPUs), digital signal processors (DSPs), etc. Lo(x, ωo, λ, t) = Le(x, ωo, λ, t) + The main advantage is that the programmer writes the ker- nel only once and OpenCL recompiles it for the compute outgoing luminance emitted luminance devices available in the system at runtime. (2) BRDF(x, ωi, ωo, λ, t) · Li(x, ωi, λ, t) · −ωi, n dωi Ω contributed luminance from all possible directions Private Private Private Private memory memory memory memory Outgoing luminance Lo, directed outward from point x in di- Work item Work item Work item Work item rection ωo in time t at wavelength λ, is the sum of luminance emitted by the object itself (Le) and the contribution of lu- minance from all possible incident directions ωi (as shown Local memory Local memory in Figure 3). Workgroup Workgroup Global + constant memory ωi1 ωo Compute device ω ωi i 3 2 Li Host memory 1 Li Figure 4: OpenCL memory model. 3 L L i o 2 Figure 4 shows OpenCL’s memory model, which is loosely Ω based on actual hardware implementations from AMD and x NVIDIA. Access to global memory is slow and should be minimized for optimal performance. Figure 3: Graphical representation of Equation 2. OpenCL’s dialect of C includes memory region qualifiers for The luminance contribution is a sum of all luminances (Li) pointers, which determine where in the memory hierarchy contributed by incident rays of light coming from the area the data they point to resides (__ of the hemisphere Ω (which extends from point x on the global, __local, __constant, __ object in the direction of the normal n). Each contribution private). It also offers a variety of specialised fixed-length vector types that map directly to SIMD registers on the un-is attenuated based on the incident angle and the BRDF derlying hardware, thus making optimisation easier. While (see subsection 1.1) of the object’s material. the language is similar to C99, it omits function pointers, bit fields, variable-length arrays, and most notably—recursion. The global illumination equation (2) has no general analyti-Instead of the standard C library, a set of math-centered cal solution and has to be solved numerically using a Monte functions is provided. Carlo approach. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 26 2. IMPLEMENTATION Random number generation is an important part of every Monte Carlo simulation. Marsaglia’s xorshift pseudorandom number generator (PRNG) was chosen for this program, as it is very fast and simple to implement on both CPU and GPU architectures [1]. Figure 6 shows the implementation of the aforementioned pseudorandom number generator. static float gen_random(unsigned int *seed) { unsigned int y = *seed; (a) 38 objects, CPU, ∼34 FPS. (b) 38 objects, GPU, ∼88 FPS. y ^= y << 13; y ^= y >> 7; y ^= y << 17; *seed = y; /* Convert random bits to float between 0 and 1 */ union { float f; unsigned int ui; } res; (c) 938 objects, CPU, ∼1.7 (d) 938 objects, GPU, ∼4.8 res.ui = (y & 0x007fffff) | 0x40000000; FPS. FPS. return (res.f - 2.0f) Figure 5: Comparison of rendering on CPU vs. GPU. * 0.5f; } The test implementation renders a sphereflake (a fractal Figure 6: PRNG implementation in OpenCL. structure made out of spheres) using the path tracing al- gorithm. Figure 5 shows a comparison of rendering speeds 3. CONCLUSION on a CPU (Intel Core i7-4771: 4 cores, 8 threads) and GPU A simple interactive path tracer capable of rendering scenes (integrated Intel HD4600: 280 compute units). The GPU is consisting of spheres was implemented using OpenCL. It about 2–3-times faster in all tested cases, as was expected. runs on both CPU and GPU architectures, with GPUs out- Adjusting the OpenCL workgroup size might give even bet- performing CPUs, as expected. ter performance. Further work might include sharing the load between CPU The program is implemented in C and OpenCL 1.0, it con- and GPU (currently, rendering is performed on either the sists of 783 lines of code in total (545 lines for the C driver CPU or the GPU, but not on both at the same time), which and 238 lines for the OpenCL kernel). The user interface should increase performance. A system for automatically is implemented with GLUT and OpenGL (OpenGL is used tuning the amount of samples rendered on the CPU vs. GPU only for drawing text and pixels into the window, all com- based on their speed would also be a useful addition. An- putation is done with OpenCL). other crucial improvement would be to use a spatial indexing structure to accelerate ray–scene intersections (currently, a Rendering is parallelised by pixels (each pixel represents naive linear traversal is used). an independent work unit). The objects, emitters (lights), and camera parameters are saved into the constant memory of the compute device, while the raw framebuffer, random 4. REFERENCES number generator state (seed) for each pixel, and the pixels [1] G. Marsaglia. Xorshift RNGs. Journal of Statistical themselves are stored in the device’s global memory. Software, 8(14), 2003. http://www.jstatsoft.org/v08/i14/paper. While the natural implementation of path tracing is recur- [2] M. Pharr and G. Humphreys. Physically Based sive, it is also possible to implement it iteratively, which is Rendering, Second Edition: From Theory To what we did. The main reason for this was because OpenCL Implementation. Morgan Kaufmann Publishers Inc., doesn’t support recursion, as it is expensive and slows down 2010. execution. [3] P. Shirley and R. K. Morley. Realistic Ray Tracing, 2nd edition. A. K. Peters, Ltd., 2003. Normally, the bouncing of rays is terminated using Russian [4] J. E. Stone, D. Gohara, and G. Shi. OpenCL: A roulette [2], however, that is not the most appropriate so-parallel programming standard for heterogeneous lution to use on GPUs. Russian roulette terminates the computing systems. Computing in science & path based on a randomly generated number, which gives engineering, 12(3):66, 2010. an uneven workload distribution. GPUs perform best when [5] C. Wynn. An introduction to BRDF-based lighting. all threads within a work unit execute the same code. If http://www.cs.ucla.edu/~zhu/tutorial/An_ a single thread has to perform a branch, it slows the other Introduction_to_BRDF-Based_Lighting.pdf. threads down. To solve this problem, a constant number of bounces is used instead. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 27 StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 28 A Garlic Clove Direction Detection based on Pixel Counting Pavel Fičur University of Primorska, UP FAMNIT Glagoljaška 8, SI6000 Koper Slovenia pavel.ficur@upr.si ABSTRACT detection is needed, in order to make sure that the garlic In agricultural communities it is well known that planting clove is inserted in correct direction into a planting mecha- garlic is very labor intensive work. The main problem con- nism. Some research was recently done in the field of com- sists of proper orientation of garlic clove at the time of in- puter vision in agriculture [2, 3]. Most research was done serting it in the soil. Because of the garlic’s irregular shape, on recognition of fruit, not on seeds. In a recent paper [1], this is usually done efficiently by human hand. As develop- authors reported a very high rate, 97%, of correctly oriented ment of technology makes use of small cameras very accessi- garlic cloves using computer vision. The method proposed ble for a wide range of applications, garlic clove orientation in the article is to collect information of root or bud [1] and detection could be efficiently implemented in garlic planting to calculate edge flexion for determining if it is the root or with machinery. A detection algorithm using simple pixel the bud. The aim of this work is to use the entire shape and counting method and an experimental result are presented. compare the information from root and bud using a sim- Experimental rate of proper direction detection was 100%. ple method of pixel counting. Authors [1] also noticed the problem of outdoor light which is solved in our case. Exper- Keywords imental work was conducted in laboratory and showed that the 100% rate of correctly detected orientation is possible. garlic planting, computer vision, direction detecting With a precise design, also on the field such result can be obtained because of the artificial light, which is crucial for the correct image. 5000 2. METHODS 4000 In this section data collection, image manipulation and de- tection algorithm are presented in detail. 3000 2.1 Collecting data 2000 For further data analysis of each clove, a black and white (bi-1000 nary) picture is needed. The picture is first taken in RGB mode and consequently set to binary (0/1) mode. In order 0 to avoid the influence of daylight, a small dark box is used. 0 50 100 150 200 250 Under the box there is an artificial light source emitting the light through an opaque glass. On this glass there is the garlic clove. The box is designed in such a way that the gar- Figure 1: Histogram of the initial capture. lic clove can enter only in right (root down) or wrong (root up) direction [1]. In this way the RGB picture is of high 1. INTRODUCTION contrast composed of near white pixels for the background and near black pixels for the clove shape. The histogram of In the process of planting garlic, proper orientation of garlic ligt level on x-axis in Figure 1 shows the gap between two cloves is crucial. The garlic cloves have an irregular shape main groups of pixels. The picture is further transformed and are of different length and width. These properties make into a matrix with 0,1 entries, 0 for background (white), 1 automated planting more difficult. An efficient orientation for foreground. The 1 are determined by the pixels covered by the clove. In Figure 2 we can se the original and bina- rized picture of the same garlic clove. 2.2 Extracting relevant elements As the garlic cloves vary in size, there may be some zero columns in the matrix, belonging to the left and the right parts of the matrix, and also there may be some zero rows, belonging to the upper and the bottom parts of the matrix. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 29 These pixels belong to background and are not relevant for our direction detection, moreover, this pixels could disturb the direction detection algorithm. All these zero columns and rows must be deleted. At this point there is a matrix representing the shape of the garlic clove without any empty (only 0) row or column. 2.3 The orientation detection algorithm As the picture is a binary matrix, it is easy to count ones in each row and adding the count to the sum of previous rows. This procedure is done twice, first from top to bottom Figure 2: Original and binarized picture of the same and second from bottom to top. If the root of the clove is clove. below, the bottom sum increases faster than the top sum and opposite for the situation where the root is up. At this moment we can select the sum results of two central rows of the clove from top and bottom, subtract the bottom sum from the top sum and check the sign. If the matix contain 4 x 10 odd number of rows, we omit the central one. If the sign is 8 negative, then the clove is correctly oriented (root below), 7 if the sign is positive, then the clove is in wrong direction 6 (root above) and must be rotated before going further in the 5 machinery. 4 3 3. EXPERIMENTAL WORK 2 All steps conducted in laboratory are explained here. At the 1 end of section the results are presented. 0 −1 3.1 Getting and preprocessing pictures −2 Three hundred garlic seed cloves were selected for testing 0 100 200 300 400 500 and inserted one by one into the box of dimensions 3 by 6 centimeters. From downside the clove was artificially illumi- Figure 3: Diagram of a correctly oriented clove. nated through the opaque glass. The seeds were of different The x-axis represents the row distance from bottom height and width and freely oriented (up, down, different or top, the y-axis represents the cumulative sum of angle). Then the RGB picture was acquired, cropped, con- both countings. The line starting and ending at zero verted in black and white mode and saved in the size of 400 level is the difference. by 800 pixels. See Figure 2. For this purpose the everyday usual photographic camera Canon G10 mounted on a tripod was used. GIMP, GNU immage manipulation program was used for image manipulation. 3.2 Direction detecting 4 x 10 10 For matrix manipulation Matlab was used. Each saved im- age was imported, counting of pixels in each row was per- 8 formed from top to bottom and from bottom to top. In Fig- ure 3 the line of the difference (bottom line) is non-positive 6 sign, that is, the clove is correctly oriented. In Figure 4 the line of the difference (bottom line) is non-negative sign, 4 that is, the clove is incorrectly oriented. Both figures also show the cumulative sums generated by algorithm. The sum 2 that increases faster is the sum of the root which is wider than the bud. For the purpose of answer, the two middle 0 rows were selected and the sign of difference of the sums was checked. For this part of work the maximum time used was −20 100 200 300 400 500 600 0,3s in Matlab on an i5 core personal computer. 3.3 Results Figure 4: Diagram of an incorrectly oriented clove. Orientation of all the three hundred cloves was correctly The x-axis represents the row distance from bottom detected, so 100% of correctly orientation detected cloves or top, the y-axis represents the cumulative sum of was achieved. The reason for success of the method is quite both countings. The line starting and ending at zero obvious because of the natural shape of garlic cloves. level is the difference. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 30 4. CONCLUSIONS In an experimental way we can conclude that the problem of detecting garlic clove direction can be solved. Now we have to pass on the field. A small camera has to be connected to some small computer, e.g. Raspberry PI and mounted on a real planting machinery. The most important data we need is the speed of planting. Is the couple camera-computer enough fast? Another question is the minimum amount of pixels for efficient run on the real machine in the sense of correct detection. 5. REFERENCES [1] G. Chi and G. Hui. Direction identification system of garlic clove based on machine vision. Telkomnika, 11(5):2323–2329, 2013. [2] Y. Wang, L. Xu, X. Zhao, and X. Hou. Maize seed embryo and position inspection based on image processing. In D. Li and Y. Chen, editors, Computer and Computing Technologies in Agriculture VII, volume 420 of IFIP Advances in Information and Communication Technology. Springer Berlin Heidelberg, 2014. [3] S. Yamamoto, K. Kobayashi, and Y. Kohno. Evaluation of a stawberry-harvesting robot in a field test. Biosystems Engineering, 15(2):160–171, 2010. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 31 StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 32 Simple approach to finding repeated patterns in opponents Texas Hold’em no-limit game Gregor Vohl Janez Brest University of Maribor University of Maribor Faculty of Electrical Engineering and Computer Faculty of Electrical Engineering and Computer Science Science Smetanova ulica 17 Smetanova ulica 17 Maribor, Slovenia Maribor, Slovenia gregor.vohl@gmail.com janez.brest@um.si ABSTRACT notes on paper. Players can make notes only in their heads. Collecting information about opponents in poker games is Finding the weaknesses of opponents is always a hard task. a hard task. You need to pay close attention to everything Opponents are always trying to improve their games and they are doing. Watch them closely and try to figure out limit their signs of weaknesses. A lot of poker players wear if they are doing anything that is abnormal - playing with sun glasses, are totally frozen, are not talking, are trying chips, breathing, shaky hands, swallowing saliva, etc. Any- to show no emotions, etc. The famous poker quote about one who has played poker before knows that this is the differ- finding the weakest player at the poker table says: ”If after ence between good and great players. It is a hard task to do ten minutes at the poker table you do not know who the at a table with living opponents but it is even harder if you patsy is − you are the patsy” [8]. Try to imagine now how do not see your opponents and can not get any information hard it would be to find the weakness of an opponent if you about the players behavior or their emotions. Poker-bots see can not see him acting behind the table and he will have no only actions from opponents and based on that they must emotional downfall. build an opponent profile of play actions. Analyzing old poker games was the key concept at imple- This article shows how it is possible to find repeated pat- mentation of poker-bot Rembrant. We used the advantage of terns in opponents games. We will show an easy approach the communication protocol from the University of Alberta. to search for important patterns in history games. With They developed an easy communication protocol which is the help of additional information a poker-bot can become used at ACPC competitions [1]. With the help of commu- a better and smarter player and win more chips from oppo- nication protocol it is very easy to have a great overview of nents. historical games. We have written a simple wrapper around the history games and implemented a searching algorithm 1. INTRODUCTION for patterns of similar play actions of opponents. Our ap- proach was to search in the history games for repeated ac- Texas Hold’em belongs to games with imperfect information. tions. Because the opponents are poker-bots it can happen That means that when a player has to make a decision some that they will repeat the same actions even if they are losing information is not available. A single player does not see the or winning more often than living players. A poker-bot can cards of other players, does not know the next card from the extend the number of chips from an opponent when a play deck and can only hope for a better end combination than an pattern is detected and repeated more frequently. opponent. This is a game of experience and intuition. The important part for the better play is analysis of old games The University of Alberta developed a poker-bot called Poki and learning from good and lost games. The key concept is [4] which uses opponent modeling with the help of a weight to find if there are repeated patterns in opponents games. table. A weight table is a group of variables collected in When playing online poker players can make notes about the a single table. Values of the table are called weights. The opponents which are available later when the same opponent weights in the table are changed according to the method is back at the table. Because there are a lot of different of modeling. The modeling technique focuses on all possi- players, notes can be useful for identifying and remembering ble starting hands because each opponent plays a different the weaknesses in opponents games. At a live poker table strategy. it is a bit complicated because it is not allowed to write The University in New Zealand developed poker modeling with the help of machine learning techniques. They used decision trees and artificial neural networks [6]. This article is organized as follows. Section 2 shows the ba- sics about the Texas Hold’em poker game. What are the poker actions and what are the end combinations which de- cide the winner. Section 3 presents our algorithm for find- StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 33 ing patterns with the help of communication protocol. A version of an algorithm with and without finding patterns Table 1: Preflop action patterns approaches is presented in Section 4. Section 5 concludes Patterns Poker actions this paper. fold to bet rf re-raise rr 2. TEXAS HOLD’EM BASICS call bet from 1st position rc limp raise cr Poker is part of card games where the end combinations are call to re-raise bet from 1st rrc ranked. Poker games are unique because of betting rounds. fold to re-raise bet from 1st rrf The game has two main parts: preflop and postflop [2]. Preflop is part of the game with pocket cards. Each player • gets two closed cards hidden to other players. There is only RE-RAISE − player pushes higher amount of chips one round of betting. From a deck of 52 cards it is possible into play than the player before him. to form 1326 different starting hands [5]. The best starting • CALL − player equals the amount of chips from the hand is pocket aces and the worst is 7 and 2 offsuit (different bet. colors). There are two manditory bets called small blind and big blind. The player next to dealer must pay the small blind • FOLD −player is out of the current game. bet and the player next to him must pay the big blind bet. Blinds are set before the start of current game. Big blind is • ALL-IN − player pushes all of the chips into the game. normally two times bigger then the small blind value. 3. COLLECTING INFORMATION ABOUT Postflop is part of a game where the community cards are OPPONENTS visible to all the players. There are 3 levels (flop, turn and river) and 3 rounds of betting. 3 community cards are added The most important fact when collecting information about on the flop and one each for turn and river. The player with opponents is how often they fold or how often they call the best combination of 5 cards (pocket cards and commu- a raise of our poker-bot. The algorithm searches through nity cards) is the winner. the history of games for patterns where the opponents were forced to fold because of bad cards or because of too high End combinations: bets and for situations where opponents called raises when a poker-bot had weak starting hands. When a pattern was discovered the poker-bot was trying to play the same com- • royal flush: A K Q J T suited, bination more often. We defined variables for some possi- ble actions and then incremented the values when they ap- • straight flush: five suited consecutive cards peared. We calculated the percentage values based on the (e.g. 4 5 6 7 8), number of played games. With the help of percentage val- • ues we defined poker-bot actions. When an opponent folds quads: four cards with the same value more frequently it is better to check/call the actions to gain (e.g. 7 7 7 7 *), more chips. If the opponent is an aggressive player then a • full house: trips plus pair good strategy is to play passive and tight and wait for good (e.g. A A A Q Q), starting hands because the amount of gaining chips will be raised. • flush: five cards of the same color, Collecting information was separated for preflop and post- • straight: five offsuited consecutive cards flop. For preflop games the key fact was how often the op- (e.g. 5 6 7 8 9), ponent called or re-raised bets from the 1st position and for • trips: three card of the same value the postflop games how often the opponent folded to bets or (e.g. 8 8 8 * *), re-raised them. • two pairs: 3.1 Collecting information on the preflop game (e.g. A A K K *), During the preflop phase of poker game the algorithm searches • pair: two cards of the same value for action patterns as shown in Table 1. First two lines of (e.g. 9 9 * * *) and Table 1 are common patterns for preflop and are possible in every single game. The next two lines are possible only • high card: (e.g. A * * * *). when poker-bot starts the preflop game. The last two lines cover three actions and the opponent must start the preflop game. Poker actions: The poker-bot Rembrant tries to play almost every preflop • CHECK − player is not betting and the right to play hand. The exception is when the starting cards are not very moves to the next player. good. In this case the poker-bot will try to minimize the amount of chips so that the flop community cards will be • BET or RAISE − player pushes desirable amount of as cheap as possible. The important fact is the percentage chips into play. of games where opponent is betting after limping from the StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 34 If the re-raise amount of opponent was bigger than 3 times Table 2: Postflop action patterns the raise of poker-bot then the poker-bot continues to play Patterns Poker actions only with strong starting hands. Calling a big re-raise with fold to bet rf weak hand means loss of chips. If the poker-bot has good fold to re-raise rrf starting cards (pocket pairs, AK or AQs) it will re-raise to call bet rc all-in. In the case of a normal re-raise, the weak starting cards are folded anyway. 1st position. If the percentage is high then the poker-bot 3.3.2 Postflop pot size Rembrant will not limp from preflop and will fold instead. A high percentage of opponent preflop raise is optimal when The raise amount in postflop is the same for all three levels starting cards are good. In this case the poker-bot wants as of postflop. The raise amount is no longer measured on the many chips as possible in the pot. raise of the poker-bot but on the pot size. Pot size is the current amount of all the chips in the play. 3.2 Collecting information on the postflop game We defined three different raise types: During the postflop phase of a poker game our algorithm searches for action patterns shown in Table 2. The patterns are equally balanced for all three postflop levels. We decided • under pot raise. to collect only the more common postflop actions here. In Table 2 are only three actions. The first two lines of the table • pot raise. present the patterns where the opponent folds to bet or re- • over pot raise. raise. The last line of the table shows the pattern where the opponent is calling the bets. This pattern can occur more than once during a single game. Under pot raise is a raise from opponent where the amount of opponents’ raise is under the current pot size. Pot size Every time before making an action in the postflop game the raise equals the amount of opponents’ raise. Over pot raise variables for the patterns are checked. The most important is a raise from opponent where the amount of opponents’ one is the call bet variable. When the percentage of calls raise is over the current pot size. is high then the poker-bot will try not to raise in case of semi-bluff situations. Instead he will try to get the next As normal raise we defined only the pot raise. That means community card for free or for low price. Aggressive play is that the opponent is raising for as many chips as in the suitable only when the fold variables have high percentages. current pot. The problems is when the opponent is betting If the fold and call variables have low percentage value then under or over the current pot size. When betting over the poker-bot will try to play tight aggressive - raise only with pot size it can mean that he is bluffing and wants us to good cards and fold in the case of weak end combination, to fold the current hand. This is very common action on the any bet of the opponent. river when acting first - it can mean monster hand or air. Variables for action patterns play an important role in this 3.3 Collecting other information scenario. If it is possible to detect the purpose of opponent’s Another very important thing to focus on is the amount of a raise with variables poker-bot has an easy choice. Other way raise by the opponent. Our algorithm calculates the average around it is very hard to predict and in most cases the better value of opponents raises through the preflop and postflop options is to fold if poker-bot has no good end combination. games. In the case of under pot raise an opponent can have a very 3.3.1 Preflop opponents average raise strong end combination and wants some more chips. To In preflop we decided to track the raise from the 1st position make under pot raise while bluffing is not a common ac-and the re-raise amount of the opponent. The average raise tion because small pot bets are normally called every time. was valid only after a certain number of games had been Poker-bot will call all raises of opponents when the under played. This is important because the average raise value pot size conditions are fulfilled. Poker-bot is calling or re- must become stable and the new values can not change the raising raises only with good end combinations. Hands with average too much. Preflop raise from the 1st position was multiple outs to hit a better end combination are called only an average raise if the following condition was met. when a normal raise from opponent is played. oppRaise−2∗BigBlind < oppRaise < oppRaise+2∗BigBlind 3.3.3 Opponent moves all-in (1) A very important fact to know is also how often an opponent If the raise amount was not within that range the poker-bot pushes all-in. All-in is a special bet or re-raise when the starts to play careful − can mean that the opponent has player moves the entire stack into the game. Losing all- strong starting hands. When the average raise was over the in means game over. Players usually play all-in with the limit the weak starting hands are automatically folded. All best cards or good end combinations. All-in move is also preflop raises from an opponent are called when the lower a good bluff method and can be performed in the preflop average raise limit is not reached. and postflop phases of the game. When putting all chips in play in the preflop game this is normally a coin flip situation Normal preflop re-raise from an opponent is defined as: and a little luck is needed to win it. We defined a variable oppReRaise > botRaise ∗ 3 (2) which is counting the number of opponents all-in moves. If StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 35 crucial. Because opponents got better and they improved Table 3: Results of Rembrant 3.0 and Rembrant 1.0 their games and it got harder to get chips from them. To vs Lucky7 protect the gained chips is therefore the key element of the Experiment Rembrant 3.0 [chips] Rembrant 1.0 [chips] poker-bots and that is why reducing the losing chips is so 1 38114 -2489 important. 2 25198 8055 3 28961 7115 5. CONCLUSION We have introduced a simple approach for modeling op- ponents actions by searching for patterns in history game. the percentage value is high that means that the opponent Variables for each pattern were converted to nominal val- is bluffing most of the time. In this case poker-bot will try ues of played games. The pattern approach is used before to call all-in moves from opponents more often with lower making the final decision about the current action the pat- ranked end combinations. tern approach is used. When a pattern is discovered current action can change in order to improve the amount of chips 4. TESTING AND RESULTS in play. In the case of losing hands poker-bot can decrease For testing we used the poker-bot Rembrant version 1.0 [7] the number of losing chips and in the opposite case it can which does not use the algorithm for finding patterns. Our improve the number of winning chips. The results from Ta- latest version of the poker-bot Rembrant is 3.0 and it uses ble 3 shows that the algorithm for patterns is effective and an approach for finding repeated patterns. We wanted to improved the previous version of poker-bot Rembrant. The test and to see whether the introduced approach improved latest version won more chips from opponents than the pre- our poker-bot. Comparison between two different versions vious one. was not enough so we also tested our algorithm against a poker-bot called Lucky7 [3] which was the result of another For future work we will try to add more pattern variables to team from our University. Lucky7 was one of the top 3 the algorithm. An increased number of variables would give poker-bots at ACPC 2011 competition. Lucky7 is a poker- a much better overview over the actions from the opponents. bot which uses multi poker-bots to define the game action. We will also try to improve the current game actions when Every poker-bot has a different style of play and a voting a play pattern is discovered to gain even more chips from system defines the current game action. opponent or to lose less chips. One of our goals is also to create a virtual simulation of the current game with all pos- The testing model was divided into two parts with 5000 sible future actions to see how a single action of poker-bot hands. After the end of the first part the poker-bots switched can effect opponents games. positions and the starting cards. All the collected informa- tion from the bots were deleted. Switch players was impor- 6. REFERENCES tant to erase the factor of luck. Small blind is 50 and big [1] Annual computer poker competition. blind is 100. Both players have a starting stack of 20.000 http://www.computerpokercompetition.org/. chips to leave enough space for combinations. Starting stack Accessed: 2014-08-10. was reset after every single hand was played. Decks of cards [2] S. Braids. The intelligent guide to Texas Hold’em poker. are generated with a random number generator. The gener- Intelligent Games Publishing, 2003. ator used the same starting seed every experiment to ensure the same deck of cards in every game. [3] B. Butolen, S. Zemlic, M. Cof, and M. Zorman. Lucky7 poker playing multi agent system. 4.1 Results http://www.computerpokercompetition.org/ downloads/competitions/2011/presentations/Team_ Table 3 shows the results of the experiment against the Lucky7_Poster.pdf. Accessed: 2014-08-29. Lucky7 poker-bot. Rembrant 1.0 was better only in the [4] A. Davidson, D. Billing, J. Shaeffer, and D. Szafron. first experiment of the test and lost the other two exper- Improved opponent modeling in poker. Proceedings of iments. All in all the Lucky7 poker-bot was better. The The 2000 International Conference on Artificial results show that the Rembrant 3.0 poker-bot improved a Intelligence (ICAI’2000), Las Vegas, Nevada, pages lot with the help of the approach for finding patterns in 1467–1473, 2000. opponents games. Rembrant 3.0 was better in all three ex- [5] R. D. Harroch and L. Krieger. Poker for dummies. periments and won a lot more chips from the opponent than Wiley Publishing, 1997. the previous version. The pattern approach helped the new [6] McNally, Patrick, Raffi, and Zafar. Opponent modeling version of poker-bot to save chips in the cases of weak hands in poker using machine learning techniques, and to gain more chips from opponents when its hands were northwestern university. good. This was possible because Lucky7 is still a case based poker-bot and its actions started to repeat after a certain http://www.cs.northwestern.edu/~ddowney/ amount of time. The algorithm detected the repeated pat- courses/349_Fall2008/projects/poker/Poker.pdf. Accessed: 2014-08-20. terns and tried to take advantage of them. With the help of the algorithm it reduced the loss and improved the gained [7] G. Vohl, B. Boskovic, and J. Brest. A rembrant poker chips by almost 4 times, to the version Rembrant 1.0. bot program. Elektrotehniski vestnik, 79(1-2):13–18, January 2012. In tests against better opponents poker-bot will not win all [8] M. J. Whitman and M. Shubik. The Aggressive the games but will reduce the amount of chips when playing Conservative Investor. Random House, New York, 1979. a losing hand. Reducing the chips when losing a hand is StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 36 Histogram of Oriented Gradients Parameter Optimization for Facial Expression Recognition Uroš Mlakar University of Maribor Smetanova 17 Maribor, Slovenia uros.mlakar@um.si ABSTRACT prise’, ’Fear’, ’Happiness’, and ’Sadness’. Each of the six This paper proposes a method for optimal parameters re- emotional states can be seen in Figure 1. garding the Histogram of Oriented Gradients (HOG) filter. Differential evolution (DE) was used for the HOG param- eter tuning using images from the MMI Facial expression database. The database was split into three parts: training, development and testing set. HOG parameters were tuned on the development set. The obtained result show an im- provement of the evolved parameters during the evolution, while the experiments on the testing set failed to provide a better performance in the end. Keywords Facial expression recognition, histogram of oriented gradi- ents, evolutionary algorithm, feature selection, differential evolution Figure 1: Six prototypic emotions [1] 1. INTRODUCTION Facial expression recognition applications can be divided Facial expression recognition has been a very interesting into two areas, depending on how the facial features are topic since the early 90s. There have been many advances extracted: geometric feature-based and appearance-based over the past decade in the areas of face detection, face methods. The main difference between the two, without go- tracking, feature selection and facial expression recognition. ing into detail, is the way information is extracted from the Facial expressions are a very powerful tool for humans to face image. With appearance based methods, a filter is usu- communicate, show emotions and their intentions, and thus ally applied to the face image, while geometric feature-based makes facial expression recognition an interesting research methods rely heavily on the shapes of important facial fea- topic in the area of human-computer interaction. Over re- tures i. e. eyes, mouth, nose, and wrinkles. In this paper cent years facial expression recognition has gained attention appearance-based methods were used for facial expression because of its wide of applications, which reach many differ- recognition. A facial expression recognition system usually ent areas like computer graphics, robotics, security, driver consists of three parts: namely face detection/tracking, fa- safety, etc. The main goal of researchers in the field is to cial feature selection and facial feature classification. In the build rich and responsive graphical user interfaces. The most first part a face is detected using the Viola-Jones detector or promising application is the application of driver attention the active appearance models. The next step after localizing control. When the driver looses attention (sleepiness) the the human face is to extract facial features. Such features are car automatically slows down [11] [6]. Regardless of the used to describe the facial expression and, consequently, the advancement of technology, facial expression recognition re- emotion of monitored person in the digital image or video mains a difficult task to this day. Researches in the area clip. The focus of this paper was the facial feature selection of facial expression recognition usually focus on recogniz-part, where an evolutionary algorithm (EA) was used to se- ing six prototypic emotional states: ’Disgust’, ’Anger’, ’Sur- lect the optimal parameters for the appearance descriptor Histogram of Oriented Gradients (HOG), to achieve maxi- mum facial expression recognition efficiency. The last part of a facial expression recognition system is a classification module, which classifies extracted facial features into one prototypic emotion (usually into six or seven emotions). The rest of the paper is organized as follows. In Section 2 some of related work is presented in the area of feature selection. In Section 3, an EA will be presented, which was used for this study, then in section 4 HOG appearance descrip- StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 37 tor will be described. In Section 5 we present the proposed As can be seen from Eq. 2, the crossover rate Cr is defined method and in Section 6 the results will be presented. With at the interval [0, 1] and it defines the probability of creating Section 7 we will conclude this paper with some findings. the trial vector parameters ti,j . The jrand index is responsible for the trial vector to contain at least one value from the 2. RELATED WORK mutant vector. After crossover, some values of the vector A number of works in parameter optimization/tuning using may be out of bounds, meaning that these values must be evolutionary algorithms have been made. In [10] a genetic mapped to the defined search space. algorithm (GA) was used for optimizing the parameters of a set of Gabor filters for the purpose of on-road vehicle detec- The next step in the evolutionary process is the selection of tion. Liu et. al [7] used GA for face recogition. The origi- the fittest individuals. During the selection process the trial nal data was projected to a lower dimension using principal vector ui, competes with vector xi from the population. The component analysis (PCA), where the rotations of the basis one with the better fitness value survives and is transfered vectors were optimized using a fitness function defined by to the next generation. performance accuracy and class separation. Yen et. al [12] used GA for optimal ellipse parameters estimation for facial 4. HOG DESCRIPTOR feature extraction. Histograms of oriented gradients (HOG) were developed by Dalal and Triggs [4] in 2005. They are based on orientations 3. EVOLUTIONARY ALGORITHM of image gradients and have been widely used for pedes- An evolutionary algorithm is a generic population-based op- trian detection. Some researchers have tried it for facial timization algorithm. It uses mechanisms which are inspired expression recognition problems and obtained very good re- by biological evolution, such as reproduction, mutation, re- sults. Orrite et al. [8] proposed a facial expression recog- combination, and selection. The solutions for the optimiza- nition system based on class specific edge distribution and tion problems are the individuals in the population, and the log-likelihood maps, the latter were used for building a list fitness function provides a metric for measuring the quality of HOG features supplemented by SVM classifier. Exper- of the solutions. In this paper an EA called differential evo- iments on the CK database pointed out 83.3% recognition lution (DE) has been used for optimal parameter selection rate. Gritti et al. [5] studied HOGs, local binary patterns of the HOG filter. Hereinafter the DE algorithm will be (LBP) and local ternary patterns (LTP) for emotion recog- described in detail. nition. Recognition rate of 92.9% on the CK database was obtained by using the LBP appearance descriptors. The idea 3.1 Differential Evolution of the descriptor is very simple, yet it provides a powerful description of the object within the image. The feature vector Differential evolution (DE) [9] is a population based opti- obtained with the HOG descriptor is calculated on a grid mization algorithm used for global optimization. It is sim- of uniformly spaced cells, which are grouped into overlap- ple, yet very effective in solving various real life problems. ping blocks, which account for contrast normalization. Let The idea of DE is a simple mathematical model, which is us take a closer look at the parameters with the following based on vector differences. example: The population of the DE algorithm consists of N p individ- Bins Cell size Block size Orientation Clip value uals xi, where i = 1, 2, ..., N p, and each vector consists of 9 16 2 0 0.2 D-dimensional floating-point encoded values xi = {xi,1, xi,2, xi,1, ..., xi,D}. In the evolutionary process, individuals are In the provided example an image is divided into cells which evolved using crossover, mutation and selection operators, are 16×16 pixels in size. Within each cell gradients are com- which are controled by a scale factor F and crossover rate puted using the filter [-1 0 1] in both directions. The orienta-Cr. The following mutation strategy was used: tions are then placed into 9 evenly distributed bins. Whether oriented gradients are used, the bins are distributed from 0◦ to 360◦ or 0◦ to 180◦ otherwise. These cells are then grouped m together into 2 × 2 blocks. In the final stage all values in i = xr1 + F ∗ (xr2 − xr3). (1) the feature vector are clipped using the clipping value. In Figure 2 HOG feaure vector calculation is presented. For For the creation of a mutant vector mi, three random vec- more information refer to [4]. tors from the population are selected, defined by indexes r1, r2 and r3. indexes are mutually different and differ- ent from i. They are selected uniformly within the range 5. PROPOSED METHOD {1, 2, ..., N p}. The scale factor F is defined within the range In this section we will briefly describe our algorithm for opti- [0, 2]. The mutant vector m mal parameter selection of the HOG appearance descriptor. i is obtained by adding a scaled vector difference to a third vector. The values of the HOG filter were described in Section 4, so we will not go into detail about them here. All parameters The trial vector t were bound by lower and upper bounds as: i is then generated using the crossover rate Cr and a corresponding vector xi from the population as: • Bin size [8, 16], • Cell size [8, 16], mi,j if rand(0, 1) ≤ Cr or j = jrand, ti,j = (2) xi,j otherwise. • Block size [2, 5], StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 38 population had the form of {xi,1, xi,2, xi,1, ..., xi,D, Fi, Cri}. In each generation right before the mutation operation, the new values of F and Cr are computed using the following equations: Fl + rand(0, 1) ∗ Fu if rand(0, 1) ≤ τ1, Fi = (4) Fi otherwise rand(0, 1) if rand(0, 1) ≤ τ2, Cri = (5) Cri otherwise. τ1 = 0.1 and τ2 = 0.1 in Eq. 4 and 5 are probabilities for adapting the values of F and Cr. Fl and Fu define the interval on which the factor value F is defined. Figure 2: HOG feature vector calculation [4] 6. RESULTS • Orientation [0, 1], In this section, the obtained results are presented, and also the experimental environment is described. • Cliping value [0, 1]. 6.1 Experimental environment All values, except the cliping value, are rounded to inte- Our method was validated on a publicly available dataset ger values at the time of evaluation. Parameter selection called the MMI Facial Expression database. The MMI database was done in the following matter. After each generation is a web searchable collection of video clips and images in of DE all individuals were evaluated in the following way: which the subjects display a wide variety of facial expres- Each image in the training set was transformed using the sions. Subjects differ in age, ethnicity and sex. There are a current individual xi of HOG parameters, giving a feature total of 69 subjects in the database. Some videos have been vector f . All vectors f were sorted according to the emo- annotated with emotions labels. For testing purposes only tion label of the images. With that in mind, after sorting, those subjects from the MMI database having emotions la- we ended up with six types of feature vectors, each belong- bels were selected. Thus, human face images from just 27 ing to a prototypic emotion. The next step was to evaluate volunteers were used in our evaluation procedure. Because the efficiencies of the current set of HOG parameters using data in the MMI database is in the form of video clips, those a classifier. The classifier used in this paper was support frames had to be extracted, where the displayed emotion was vector machines (SVM) using a radial basis function kernel. at the peak (apex) throughout its duration. In other words, An implementation of SVM from the LibSVM toolbox [3] just those frames where emotions were the most expressive was used. A 5-fold-cross validation was carried out in order were selected in our analysis (3 frames). to find optimal parameters for SVM and our classification problem. Six SVMs were built, each to distinguish one emo- All the acquired frames were then split into three subsets: tion from the others. All images in the development set were training, development and testing set. By splitting the database transformed to feature vectors, and then classified using the attention was paid that no subject appeared in two sets at trained SVMs. The final output of the evaluation procedure the same time (all three sets were unique). That gave us a of a single HOG parameters’ set is a metric called efficiency total of 10 subjects for the testing set, 10 for the training of EA (EEA) (Eq. 3), which is defined as: set and 6 for development. In Table 1 we can see the total number of images per emotion in the whole dataset, while Table 2 presents a total of images per splited dataset. Nd EEA = ∗ 100, (3) Ad Table 1: Number of images per emotion in the whole dataset where Nd is the number of correctly classified images and Dis. Ang. Sur. Fea. Hap. Sad. Σ Ad the total number of all images in the development set. 102 108 132 105 132 108 687 Lastly we must also mention the jDE strategy which was The split was done randomly, with 40 % of all subjects se- used in this paper. As mentioned in Section 3, the DE al- lected for training(10), another 40 % for testing(10) and 20 gorithm has three control parameters: Np, F and Cr. The % for the development set(6). values of these control parameters are usually problem de- pendent, and their values adapt during the evolutionary pro- cess. A self-adaptive control strategy was used as in [2] for 6.2 Experimental results the control of values F and Cr. With the inclusion of con- The following settings of the DE were used for the experi- trol parameters during the evolution, each individual in the ments in this paper : N p was set to 20, G to 100, while the StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 39 7. CONCLUSION Table 2: Number of images in each dataset This paper proposed a method for optimal HOG descrip- Training Development Test tor parameters selection using the self-adaptive DE. It was 249 147 291 tested on a publicly available Cohn Kanade dataset, which was separated into three sets: testing, training and devel- opment. The HOG parameter tuning was done with the F and Cr were set to 0.9 and 0.1 respectively at the start of help of the training and development sets. The results show the evolutionary process. The N p and G were set to smaller an improvement in selecting the appropriate HOG param- values so the results could be obtained in reasonable time. eters on the development set. Due to fairly large feature In each generation N p ∗ 6 classifier learning was done. vectors, the learning part of the proposed method was very time consuming, therefore smaller values for N p and G were In Figure 3 the results of the best, worst, median and mean selected. The obtained results from tuning were also tested are reported for each generation on the development set . on the testing set. A comparative study was made between It can be seen from the graph that the DE algorithm was the best and worst results of the tuning parameters in the successful in tuning the optimal parameters on the develop- first and last generations on the test set. Because of the ment set. For the final experiment we also tried the tuned spontaneous emotions displayed by subjects in the MMI Fa- HOG parameters on the test set and compared them to the cial expression database the results in the case of the best randomly generated parameters in the first generation. The tuned parameter in the last generation gave worse results results of the comparison are collated in Table 3. than those in the first generation. In the future we would like to address the problem of spon- Table 3: Comparison between the first and last gen- taneous emotions during HOG filter parameter optimiza- eration parameters on the test set 1st generation 100th generation tion, and also test the proposed algorithm on more available datasets. best 394.37 389.32 worst 100.00 340.70 8. REFERENCES After testing the tuned parameter sets on the testing set of [1] Six prototypic emotions. the MMI database the experiment showed intersting results. https://managementmania.com/en/six-basic- The best parameter set in the last generation gave worse re- emotions. Accessed: sult then the best in the first generation, in spite of giving 2014-07-18. better results on the development set. We investigated the [2] J. Brest, S. Greiner, B. Boskovic, M. Mernik, and cause of this problem and found an answer in the database V. Zumer. Self-adapting control parameters in itself. The MMI Facial expression database mostly contains differential evolution: A comparative study on videos of subjects which are performing spontaneous emo- numerical benchmark problems. Evolutionary tional responses, that we try to classify as six emotional Computation, IEEE Transactions on, 10(6):646–657, states. The problem at hand is that each subject expresses 2006. emotions in their own particular way. Because the training, [3] C.-C. Chang and C.-J. Lin. LIBSVM: A library for development and test sets were completely independent of support vector machines. ACM Transactions on each other, and because of the way people show spontaneous Intelligent Systems and Technology, 2:27:1–27:27, emotions, the obtained results are slightly worse than those 2011. Software available at on the development set. urlhttp://www.csie.ntu.edu.tw/ cjlin/libsvm. [4] N. Dalal and B. Triggs. Histograms of oriented gradients for human detection. In Cordelia Schmid, Stefano Soatto, and Carlo Tomasi, editors, International Conference on Computer Vision & Pattern Recognition, volume 2, pages 886–893, INRIA Rhône-Alpes, ZIRST-655, av. de l’Europe, Montbonnot-38334, June 2005. [5] T. Gritti, C. Shan, V. Jeanne, and R. Braspenning. Local features based facial expression recognition with face registration errors. In Automatic Face & Gesture Recognition, 2008. FG’08. 8th IEEE International Conference on, pages 1–8. IEEE, 2008. [6] Qiang Ji and Xiaojie Yang. Real-time eye, gaze, and face pose tracking for monitoring driver vigilance. Real-Time Imaging, 8(5):357–377, 2002. [7] C. Liu and H. Wechsler. Evolutionary pursuit and its application to face recognition. Pattern Analysis and Machine Intelligence, IEEE Transactions on, Figure 3: Graph of best, worst, mean and median 22(6):570–582, 2000. for each generation during tuning on development [8] C. Orrite, A. Ga˜ nán, and G. Rogez. Hog-based set decision tree for facial expression classification. In StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 40 Pattern Recognition and Image Analysis, pages 176–183. Springer, 2009. [9] R. Storn and K. Price. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. Journal of global optimization, 11(4):341–359, 1997. [10] Z. Sun, G. Bebis, and R. Miller. On-road vehicle detection using evolutionary gabor filter optimization. Intelligent Transportation Systems, IEEE Transactions on, 6(2):125–137, 2005. [11] Mohan M Trivedi, Tarak Gandhi, and Joel McCall. Looking-in and looking-out of a vehicle: Computer-vision-based enhanced vehicle safety. Intelligent Transportation Systems, IEEE Transactions on, 8(1):108–120, 2007. [12] G. G. Yen and N. Nithianandan. Facial feature extraction using genetic algorithm. In Evolutionary Computation, 2002. CEC’02. Proceedings of the 2002 Congress on, volume 2, pages 1895–1900. IEEE, 2002. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 41 StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 42 Constructing domain-specific semantic dictionaries to supplement domain-specific knowledge bases ∗ Goran Hrovat Milan Ojsteršek University of Maribor University of Maribor Faculty of Electrical Engineering and Computer Faculty of Electrical Engineering and Computer Science Science Smetanova 17, 2000 Maribor Smetanova 17, 2000 Maribor Slovenia Slovenia goran.hrovat@um.si milan.ojstersek@um.si ABSTRACT constructing BabelNet, a wide-coverage multilingual knowl- Semantic dictionaries as well as knowledge bases are im- edge resource. The method is consisted of two steps, where portant source of information for natural language process- the first step produces a mapping between a multilingual en- ing. Using corpus and algorithms for constructing semantic cyclopedic knowledge repository (Wikipedia) and a compu- space, we can quickly construct semantic dictionary, which tational lexicon of English (WordNet). The second step fur- is exploit to supplement the knowledge base. Algorithm for ther use a machine translation system to collect a very large constructing semantic space, COALS was chosen. Semantic amount of multilingual concept lexicalization. Morsey et distance between terms in semantic space reveals their sim- al. [9] constructed DBpedia which continually extracts struc- ilarity, which is used to add semantic relationships in the tured information from Wikipedia and store them in files dictionary. Semantic dictionary also serves to assist in iden- as RDF triples. This article proposes the method for con- tifying entities and in relation extraction, to supplement the structing domain-specific semantic dictionary using domain-knowledge base. specific corpus or collection of documents in order to sup- plement domain-specific knowledge base. Keywords 2. ISSUES AT CONSTRUCTING DOMAIN- semantic space, correlation, semantic dictionary, knowledge base SPECIFIC SEMANTIC DICTIONARIES The main requirements for constructing domain-specific se- mantic dictionary are domain-specific knowledge and domain- 1. INTRODUCTION specific corpus. For domain-specific corpus bachelor thesis, Semantic dictionaries play important role in tasks of natural journals etc. can be used. For domain-specific dictionary a language processing (e.g., named entity recognition [1], rela- bag of words from corpus are considered and further ana- tion extraction [2], information extraction [3] [4]) as well as lyzed. All words are morphologically tagged and stop words knowledge bases. In semantic dictionary relations between such as pronouns and conjunctions are excluded. The most words are defined, such as hypernym, hyponym, antonym, important feature of semantic dictionaries is relations be- synonym, etc. Knowledge base is even more complex, where tween words. The examination of n2 word pairs in case of any possible relation can be defined between terms to con- n words need to be done to define relations (e.g. hyper- struct graph. It is usually represented in RDF format. An nym) between all word pairs. Because of such amount of example of semantic dictionary is WordNet [5] and an ex- pairs manual examination is unfeasible and some automatic ample of knowledge base is DBpedia [6]. Some approaches or semi-automatic approach is necessary. for dictionary construction have already been proposed in the past. Xu et al. [7] constructed a medical dictionary of 3. CONSTRUCTING DOMAIN-SPECIFIC DIC- disease terms from randomized clinical trials (RCT) using an automated, unsupervised, iterative pattern learning ap- TIONARY USING SEMANTIC SPACE proach. Navigli et al. [8] developed an automatic method for To reduce the n2 problem space, algorithms for construct- ing semantic space is incorporated in our method. Semantic ∗Corresponding author space is high dimensional space where each word is placed based on its meaning. Words in this space close to each other by specific measure (e.g. Euclidean distance, cosine dis- tance, correlation distance) are therefore semantically close, which can be exploit for efficiently constructing semantic dictionary. To place words in semantic space one of the se- mantic algorithms is performed on corpus. The most known algorithms for constructing semantic space are HAL [10], LSA [11], Random Indexing [12], COALS [13] and some al- gorithms based on dictionary WordNet [14], which do not use corpus. In our case COALS is used on bachelor theses StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 43 of University of Maribor. All documents had been processed T − w w j a,j i i,b w where unimportant part of texts (e.g. table of contents, ref- a,b = wa,j (T − wa,j ) wi,b(T − wi,b) erences) were excluded. j j i i (2) 3.1 Corpus preprocessing 3. Each cell in the matrix is further transformed with Domain-specific corpus is required for using COALS (Corre- Equation 3. lated Occurrence Analogue to Lexical Semantic). In our case the bachelor theses from University of Maribor from the field 0, x < 0 ai,j = √ (3) of computer science are used. Bachelor theses are in PDF x x ≥ 0 format, loaded in digital library of University of Maribor by students. Preprocessing of the corpus was performed in next 4. Last step is dimensionality reduction of vectors which phases: is performed using SVD (singular value decomposi- tion) [18] shown in Equation 4. In our case dimension was reduced to 800. 1. Converting PDF format to plain text using Apache Tika [15]. M = U ΣV ∗ (4) 2. Segmentation of each the bachelor thesis to first page, 3.3 Semantic similarity table of content, table of figures, content and refer- Semantic similarity between words is calculated using their ences. semantic vectors. For calculating similarity between two se- 3. Tokenization of the documents’ contents. mantic vectors well known similarity measures are used. In case of HAL Euclidean distance is used, for LSA cosines dis- 4. Morphologically tagging and lemmatization of the doc- tance is used and in our case for COALS correlation distance uments’ contents using TreeTagger [16] learned on mor- is used on the two vectors presented in Equation 5. phologically tagged corpus FidaPlus [17]. (ai − ¯ a)(bi − ¯ b) S (5) 5. Only semantically meaningful words, meaningful verbs, a,b = (ai − ¯ a)2 (bi − ¯ b)2 nouns, adjectives, adverbs and abbreviations are pre- served. Other morphological types do not reflect se- mantics of the text and are therefore excluded. Equation 5 returns similarities between -1 and 1 where 0 means no correlation between words and values 1 or -1 mean 6. Converting to lower case letters. strong correlation between words. On Fig. 1 words are pre- sented in two dimensional semantic space. For representa- 3.2 Calculating semantic vectors using algo- tional purposes dimension is reduced from 800 to 2 using MDS (multidimensional scaling) [19] where correlation sim- rithm COALS ilarity is used as a distance measure. Three different groups Semantic vector of the word represents the word in the se- are marked for clear presentation of the method. Words mantic space. Semantic vector is also context vector, be- mac and widows are close together, which means that they cause it contains encoded context of the word. Algorithm appear in similar context, but not necessary together. COALS calculates semantic vectors in the next phases. 1. Construction of the co-occurrence matrix, where each row represents observed word and each column repre- sents context word frequency. This matrix is called word by word matrix in comparison with more known LSA where word by document matrix is constructed. For observed word weights of context words are incre- mented according to its occurrence in text. In our case window of four words is used, which means that next the nearest context word get weight of 4, one context word away get weight of 3, etc. In contrast, LSA uses whole document window where all words are weighted equally. Number of columns or dimensionality of se- mantic vectors depends on the number of unique words in corpus. Because the number of unique words can be large we can limit it to 100,000 most frequent words and thus limit the dimension to 100,000. 2. The matrix is then transformed using correlation nor- malization to reduce the impact of high frequencies. Correlation normalization is calculated with Equation 2. Figure 1: Visualization of the words in two dimen- T = wi,j (1) sional space i j StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 44 3.4 Preparation of the results of COALS 4.1 Named entity recognition Algorithm COALS is used to compute semantic vectors where Named entity recognition can further improve the semantic correlation similarity metric gives us similarity between the dictionary. It is useful for defining hypernyms and is also words. Values closer to 1 or -1 mean stronger similarity prerequisite for relation extraction. Named entity recogni- whereas 0 means no similarity. Using this information the tion is used for recognition of people, animals, location, time data can be prepared for constructing semantic dictionary etc. If dog is recognized as an animal, can be that inserted where experts can define type of relations for their domain. in the domain specific knowledge base. Word pairs is then ordered by correlation absolute value in descending order so an expert can consider the most similar 4.2 Relation extraction word pairs first. Construction of the semantic dictionary is Relation extraction enables us to find relations between the much faster using this method. words from the corpus, which are then inserted into the knowledge base. For relation extraction patterns are needed 3.5 Construction of the semantic dictionary and which define relations. Examples of the patterns for relation relations ”is husband of” is shown in Table 1. An expert defines relations between the most interesting word pairs in the semantic dictionary. The interestingness Table 1: Examples of the patterns for relation ”is is calculated with the algorithm COALS, what significantly husband of” used for relation extraction reduces the number of interesting word pairs. An expert Pattern Relation ”is husband of” has some definitions predefined however his own can also be A is husband of B A ”is husband of” B added. Predefined relations are: A has wife B A ”is husband of” B A is married with B if A is male then • A ”is husband of” B synonym (two words are synonyms, when they have else same meaning) B ”is husband of” A • antonym (two words are antonyms when they have op- A and B went on if A is male then posite meaning) the honeymoon A ”is husband of” B else • hypernym (word A is a hypernym of word B, if A has B ”is husband of” A more general meaning) Semantic dictionary can be used for manually defining the • hyponym (word A is a hyponym of word B, if A has patterns for relation extraction, which can serve us for im- more specific meaning) proving knowledge base. • meronym (word A is a meronym of word B when A is a part of B, e.g., CPU is a meronym of computer) 5. CONCLUSION The paper presents the method for constructing domain spe- • holonym (inverse relation of a merony, e.g., computer cific dictionaries using domain specific corpuses. The use of is a holonym of CPU) the algorithm COALS for constructing semantic dictionary is presented. The COALS is crucial for reducing the number • troponym (word A is a troponym of word B when A of word pairs, which an expert need to consider to define more specifically expresses an action of B, e.g., typing relations between them. On Fig. 1 results of the COALS is a troponym of writing) is shown for some interesting words. In subsection 3.5 some • relations are defined however an expert is encouraged to cre- coordinate word (word A and B are coordinate words ate their own. In section 4 the procedure for supplementing to each other, when they have same hypernym, e.g. knowledge base is presented using named entity recognition CD and DVD are coordinate words, because they have and relation extraction. For further research on supplement- same hyperny medium) ing knowledge base, methods of machine translation can be • related word (when expert can define any of the above incorporated. mentioned relations, then he can define his own rela- tion or set only that word A and B are related words) 6. ACKNOWLEDGEMENT We thank the authors of the FidaPLUS corpus, who grant 4. SUPPLEMENTING DOMAIN SPECIFIC us its use for morphological tagging and lemmatization. KNOWLEDGE BASE Domain specific knowledge base contains only knowledge 7. REFERENCES from the single domain. Knowledge base is represented with [1] D. Nadeau and S. Sekine, “A survey of named entity an ontology where also rules are defined used for reasoning. recognition and classification,” Lingvisticae Knowledge base can also contain relations from semantic Investigationes, vol. 30, no. 1, pp. 3–26, 2007. dictionary, which can be supplemented from the corpus us- [2] A. Schutz and P. Buitelaar, “Relext: A tool for ing the methods of natural language processing (e.g., named relation extraction from text in ontology extension,” in entity recognition, relation extraction). For each relation its The Semantic Web–ISWC 2005, pp. 593–606, frequency can be also acquired using the corpus. Springer, 2005. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 45 [3] E. Tutubalina and V. Ivanov, “Unsupervised approach to extracting problem phrases from user reviews of products,” [4] A.-M. Popescu and O. Etzioni, “Extracting product features and opinions from reviews,” in Natural language processing and text mining, pp. 9–28, Springer, 2007. [5] G. A. Miller, R. Beckwith, C. Fellbaum, D. Gross, and K. J. Miller, “Introduction to wordnet: An on-line lexical database*,” International journal of lexicography, vol. 3, no. 4, pp. 235–244, 1990. [6] “Dbpedia.” http://dbpedia.org. Accessed: 2014-09-25. [7] R. Xu, K. Supekar, A. Morgan, A. Das, and A. Garber, “Unsupervised method for automatic construction of a disease dictionary from a large free text collection,” in AMIA Annual Symposium Proceedings, vol. 2008, p. 820, American Medical Informatics Association, 2008. [8] R. Navigli and S. P. Ponzetto, “Babelnet: The automatic construction, evaluation and application of a wide-coverage multilingual semantic network,” Artificial Intelligence, vol. 193, pp. 217–250, 2012. [9] M. Morsey, J. Lehmann, S. Auer, C. Stadler, and S. Hellmann, “Dbpedia and the live extraction of structured data from wikipedia,” Program: electronic library and information systems, vol. 46, no. 2, pp. 157–181, 2012. [10] K. Lund and C. Burgess, “Producing high-dimensional semantic spaces from lexical co-occurrence,” Behavior Research Methods, Instruments, & Computers, vol. 28, no. 2, pp. 203–208, 1996. [11] T. K. Landauer, P. W. Foltz, and D. Laham, “An introduction to latent semantic analysis,” Discourse processes, vol. 25, no. 2-3, pp. 259–284, 1998. [12] M. Sahlgren, “An introduction to random indexing,” in Methods and applications of semantic indexing workshop at the 7th international conference on terminology and knowledge engineering, TKE, vol. 5, 2005. [13] D. L. Rohde, L. M. Gonnerman, and D. C. Plaut, “An improved method for deriving word meaning from lexical co-occurrence,” Cognitive Psychology, vol. 7, pp. 573–605, 2004. [14] W. D. Lewis, “Measuring conceptual distance using wordnet: the design of a metric for measuring semantic similarity,” 2001. [15] “Apache. tika..” http://tika.apache.org/. Accessed: 2014-09-25. [16] H. Schmid, “Probabilistic part-of-speech tagging using decision trees,” in Proceedings of the international conference on new methods in language processing, vol. 12, pp. 44–49, Manchester, UK, 1994. [17] “Fidaplus.” http://www.fidaplus.net. Accessed: 2012. [18] G. Golub and W. Kahan, “Calculating the singular values and pseudo-inverse of a matrix,” Journal of the Society for Industrial & Applied Mathematics, Series B: Numerical Analysis, vol. 2, no. 2, pp. 205–224, 1965. [19] J. B. Kruskal, “Nonmetric multidimensional scaling: a numerical method,” Psychometrika, vol. 29, no. 2, pp. 115–129, 1964. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 46 Am I overtraining? A novel data mining approach for avoiding overtraining Iztok Fister Jr. Goran Hrovat Samo Rauter University of Maribor University of Maribor University of Ljubljana Faculty of Electrical Faculty of Electrical Faculty of Sport Engineering and Computer Engineering and Computer Gortanova 22, 1000 Ljubljana Science Science Slovenia Smetanova 17, 2000 Maribor Smetanova 17, 2000 Maribor samo.rauter@fsp.uni-lj.si Slovenia Slovenia iztok.fister1@um.si goran.hrovat@um.si Iztok Fister University of Maribor Faculty of Electrical Engineering and Computer Science Smetanova 17, 2000 Maribor Slovenia iztok.fister@um.si ABSTRACT his/her achievement with other people. This stage naturally Overtraining is one of the biggest problems in process of leads to a competition. To compete in sport and fight for a sport training. Especially, freshmen’s and amateur athletes medals is not easy anymore. Firstly, athletes have a devotion who do not have enough knowledge about behavior of their to the sport and in line with this assume themselves tremen- body, training and who do not have personal trainers en- dous amount of work. Secondly, they have to build/create counter overtraining. So far some theoretical and practical a good and effective training sessions. To finish race or to solutions for avoiding overtraining have been arisen. In this win a race is a huge difference and this difference is usu- paper, we propose a novel automatic solution which helps ally reflected in training, eating and resting. When amateur athletes to avoid overtraining. It is based on data mining athletes start to deal with a sport, they do not have enough which is applied to athlete’s workout history. First experi- knowledge about these three components of each prepara- ments and simulations showed the promising results. tion for sport competitions: Keywords • Training: Usually, the following set questions is asked data mining, apriori, sport training, overtraining by athletes. What is a proper training? When I have to train? How much to train? What kind of training 1. INTRODUCTION should I train? Are intervals better than long distance Recently, an expansion of massive sport events and competi- training? Where should I train? With whom should I tions have been emerged. These events capture usually sport train? disciplines, like running, cycling and triathlon [1]. Conse- • Eating: The following question are to be answered in quently, more people joined to a sport and healthy lifestyle. line with this. What I have to eat? Why I need to In fact, a lot of health and medicine organizations in the eat so much carbohydrates? Do I need to increase world also encourage people to participate in such events and vegetables and fruits intake? Is it good to drink any live a healthy lifestyle. When people start to deal with sport, powders or isotonic drinks? Why do I need to avoid they do not have any special goals. There, only one goal is eating sugar? Is it good to take additional proteins presented, i.e., I have to finish a specific sport race/event into my food? Why eating a lot of fish is useful? in which I invested so much effort! However, when someone finish more races, he/she want to go forward and compare • Resting: Here, the following questions could be emerged, e.g., Do I need to take a rest every week? When I need rest and free days? What is a good resting? Can I still have a simple and short run/cycle? Swimming helps for relaxation of muscles and body? These three components define well trained athlete. Train- ing and resting are equally important. If you train a lot you get overtraining. On the other hand, you cannot be trained enough and suffer in the races. Finally, eating has a great StuCoSReC Proceedings of the 2014 1 st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 47 impact to athlete’s performance. It is important for recovery Physiological Psychological and regeneration. Nowadays, overtraining is still one of the Higher resting heart Sleep disturbances biggest problems in sport. Almost every freshman encounter rate this problem. When freshman encounter this problem, they Changes in normal Loss of self- usually suffer for psychological problems in mind. Moreover, blood pressure confidence it is very hard to overtake this problem. So far many theo- Delayed return to Drowsiness and apa- retical, practical and abstract solutions have been developed normal heart rate thy to deal with this problem. In this paper we propose a new Elevated basal Irritability solution to help athletes to avoid overtraining. Data min- metabolic rate ing methods are applied on athlete’s workouts which were Elevated body tem- Emotional and moti- created by sport trackers. Our solution send an alarm to perature vational imbalance athletes before overtraining. Therefore, they can change a Weight Excessive, prolonged training or take a more rest and prevent themselves from loss/excessive thirst weariness overtraining. The paper is structured as follows: in the Impedded respiration Lack of appetite next section we describe some basics about sport training. Subcostal aching Fatigue The third section outlines overtraining, while section 4 re- Bowel disorders Depression views current overtraining prevention methods. Section 5 Anxiety describes data mining methods. Section 6 is devoted to the Anger/hostility sport trackers. Section 7 proposes a novel solution and the Confusion last section conclude the paper with some remarks for future work. Table 2: Physiological and Psychological indicators of overtraining syndrome 2. SPORT TRAINING Definition of Sports Training is based on scientific and ped- Causes agogical principles of planning the systematic activities or 1 Length of the competitive season processes in order to achieve the highest results in the chosen 2 Monotony of training sports discipline. The final effect of the systematic training 3 Feelings of claustrophobia process can be manifested as 4 Lack of positive reinforcement 5 Feelings of helplessness 1. athlete’s well sport form, 6 Abusiveness from authorities 7 Stringent rules 2. the increased capacity of the athlete’s organism, or in the worst case 8 High levels of competitive stress 3. overtraining syndrome. Table 3: Causes of overtraining syndrome as pro- posed in [4] Sport form can be described as a phenomenon of short-term increased capacity of the athletes organism in relation to the expected competitive capacity that is perceived on the As it can be seen from Table 2, there are many indicators subjective level. In such conditions athlete feels that he/she telling the athletes are overtrained. However, better than overcomes a certain load of sport activity with a little effort curing overtraining is avoiding overtraining. In [4], authors or that he/she is able to overcome a higher load at the max- presented some causes of overtraining syndrome. imum effort [2]. In simple terms, it means ”to be in the best shape at the right time.” Efficiency of the process of sports A lot of research were done on overtraining in the past, e.g., training also means that the athlete rationalize the time de-in [5, 6, 7, 8, 9]. voted to the training [3]. All this is happening in a real process of sports training, where in order to achieve sports 2.2 Current methods of avoiding overtraining form, athletes include different methods in the process of The main herald of overtraining syndrome is a hard training sport training and intensity of workload. Also the rest pe- and an inadequacy of the rest. People transferred training or riod is very important. For instance, Table 1 presents an competition stress differently. Athlete tolerance on training example of proper training two weeks before half marathon stress varies throughout the season periods. In line with this, race. a training process should be adapted and varied through the season period. Especially, the right strategy in the compe- 2.1 Overtraining tition period is very important. Too much high intensity Overtraining is a special state of an athlete occuring when workout with too short rest period may results in bad re- someone undergoes a bigger amount of intensity training sults at the competition or in the worst case may leads to and consequently body is not able to recover. In such state, overtraining. To avoid overtraining syndrome is not so easy. athletes are not able to compete on high level because their Athletes are most of the time very tired. It is hard to distin-body is exhausted. Overtraining leads to overtraining syn- guish when the first signs of fatigue happen. Prevention re- drome. Overtraining syndrome manifests itself in physiolog- quires good nutrition, complete hydration, and rest periods ical and psychological symptons and has an effect on ath- between exercises themselves [7]. The effective and good or- lete’s performance [4]. These physiological and psychological ganize periodization process is necessary to ensure adequate indicators as proposed in [4] and presented in Table 2. adaption of the athlete organism to the requirements of the StuCoSReC Proceedings of the 2014 1 st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 48 DAY METHODS DURATION INTENSITY 1 Race 2 Easy jogging with some accelera- 30 min Low tion 3 Rest day 4 Running (interval training 2 x 10 45 min Low & High min fast; between sets 5 min - 8 min of easy jogging) 5 Rest day 6 Running (interval training 4 x 4 min 45 min Low & High fast; between sets 3 min - 4 min of easy jogging) 7 Running (interval training 2 sets of 45 min Low & High 3 - 5 x 90 sec accelaration) 8 Rest day 9 Rest day 10 Endurance running 90 min Low 11 Endurance running 75 min Low 12 Rest day 13 Endurance running 105 min Low 14 Running (interval training 2 sets of 75 min Low & High 3 x 2 min very fast with 1 min re- covery; between sets 5 min - 8 min of easy jogging) Table 1: Example of training for 21 km race (last 14 days before the competition) competition in chosen sport discipline. mining (DM) is a method, where the main goal is extrac- tion of information from a data sets and conversion of these To avoid overtraining syndrome it can help an individual data to a form understandable by humans. The more fre- monitoring of: quently used methods in data mining are: clustering, clas- sification, regression, association rule mining. Data mining methods involve also some machine learning methods like • achievements, decision trees, logistic regression, support vector machines, etc. But nowadays also some nature-inspired algorithms [12, • mood states and 13] are employed in the data mining field. In data mining, • heart rate variability. a broad spectrum of many applications has been developed, e.g., these methods have been used for weather prediction, fraud detection, sensor data mining, surveillance and many To overcome the overtraining syndrome it can help: more. • low intensity workout, rest and relaxation strategy; 3.1 Association rule mining • exercise with very short high intensity sprints with long One of the data mining tasks is association rule mining. It rest/low intensity period, the confidence to progress was discovered by Agrawal and Srikant [14, 15] who devel- well. oped Apriori. Other algorithms for association rule mining also emerged e.g. FP-growth [16], which perform faster. As- sociation rule mining is a task of discovering rules from set 2.3 Monitoring the sport activities of transactions. The well known example is basket market Recently, training sessions of athletes are monitored with analysis however association rule mining can be applied for sport trackers. Sport trackers are a combination of software other datasets such as sports data. Each transaction is con- and hardware which are used for tracking activities during sisted of items and is also called item-set. Discovered rules sport workouts. Hardware is a smart phone, while software are the form X ⇒ Y , where X is subset of items and Y are applications which run on smart phones. For a more is usually only one item. As a result many rules are dis- detailed description of sport trackers, readers are invited to covered, for this reason a lot of measures are proposed to read papers [10, 11]. evaluate and rank these rules. The most known interesting- ness measures are support and confidence [14], others are 3. DATA MINING OF MONITORED DATA lift [17], all-confidence [18], conviction, leverage, χ2 [19] etc. Data obtained during the sport session of athletes can be Support is proportion of transactions containing items of a transferred to personal computers in a form of TCX data rule and confidence is defined as conf (X ⇒ Y ) = supp(X ⇒ sets. Normally, these data sets are analyzed using differ- Y )/supp(X). In our case transactions is consisted of items ent data mining methods. In computer science, the data e.g. training distance, average heart rate, motivation, eat- StuCoSReC Proceedings of the 2014 1 st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 49 ing, sleeping, etc., and an example of rule which can be 4.4 Applying Apriori algorithm discovered is ( M OT IV AT ION , EAT IN G ) ⇒ In this step, Apriori algorithm is launched using the pre- ( SLEEP IN G , 1.000), where confidence and support are 1 pared data (transactions). Apriori algorithm discovers rules and means that everyone who well eats and is motivated also as a result. For example, a rule ( M OT IV AT ION , well sleeps. EAT IN G ) ⇒ ( SLEEP IN G , 1.000) means that athlete who has a motivation and can eat well, can also sleep well. 4. PROPOSED METHOD FOR AVOIDING THE OVERTRAINING 4.5 Interpretation of results Our proposed method consists of the following steps: The last step in the proposed approach is an interpretation of results. From the rules which we get with Apriori algorithm, we want to get a real situation. • definition and design of framework, • discretization, 5. EXPERIMENTS AND SPORT INTERPRE- • TATION OF RULES creating transactions, The first experiments of the proposed method for predict- • applying Apriori algorithm and ing the overtraining were conducted on a real dataset. This dataset was produced during the training sessions of a pro- • interpretation of results. fessional mountain biker. The total number of workouts/transactions were limited to 50. In Figure 1, some In the remainder of the paper, these steps are illustrated in transactions are presented that were used in our experi- detail. ments. 4.1 Definition and design of framework At first, a framework needs to be defined that enable an athlete to store, maintain and use data. Let us suppose that after every training session athlete uploads a file in TCX format. This file is automatically parsed and data (total duration, total distance, average heart rate, max heart rate) are stored into the database. During uploading the file, Figure 1: Example of transactions athlete needs to answer to the five predefined questions that characterize current training session. These questions are: From experiments, interesting results were obtained. For instance, the rules are inferred from the transactions, as follows. 1. Did athlete train intervals? (Possible answer values: YES, NO) 1. Rule: ( EAT IN G , ) ⇒ ( M OT IV AT ION , 1.000) In- 2. Did athlete feel any special pains during the training terpretation: From this rule we can see that proper eating session? (Values: YES, NO) is connected with motivation. If we eat a lot we also have a motivation. 3. Do athlete have a great motivation for future train- 2. Rule: ( N O M OT IV AT ION , ) ⇒ ings? (Values: YES, NO) ( N O SLEEP IN G , 1.000) If athlete does not have a mo- tivation, he can not sleep. It might means that he is over- 4. Did athlete sleep well today? (Values: YES, NO) trained. 5. Can atlete eat normally? (Values: YES, NO) 3. Rule: ( LON G RIDE , ) ⇒ ( N O IN T ERV AL , 1.000) From this rule we can see what is usually in practice. If we ride a long rides, we do not do intervals, but going the When athlete answers to these questions, he/she is linked distance only. with today’s workout. 4. Rule: ( SLEEP IN G , ) ⇒ ( M OT IV AT ION , 1.000) Here we see can again that sleeping and motivation are con- 4.2 Discretization nected. If athlete can sleep well, then he is very motivated When all data are collected a discretization process is per- and for sure not overtrained. formed. In this process, numerical values have to be dis- 5. Rule: ( N O M OT IV AT ION , ) ⇒ cretized, since Apriori algorithm works only with categorical ( N O EAT IN G , N O SLEEP IN G , 1.000) In this rule we values. An example of discretized attributes is presented in see that athlete who does not have a motivation, can not eat Table 4. and sleep. In this case, he is very overtrained. 6. Rule: ( LON G RIDE , N O EAT IN G ) ⇒ 4.3 Creating transactions ( N O IN T ERV AL , 1.000) This rule defines long ride and After attributes are discretized, a transactions are created. that athlete can not eat. We can also see that in this case In our case, we consider one transaction for one training. athlete were not doing an intervals. Therefore, after 2 months of training, there will be around 45 transactions. More transactions we have, more accurate From previous rules, it can be seen some rules which were rules can be inferred and in line with this also the overtrain-discovered from a real dataset. In some cases we can see how ing can be predicted more precisely. to determine overtraining in rules. However, after running StuCoSReC Proceedings of the 2014 1 st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 50 Attribute Current values Discretized values How to discretize? Distance Numerical (km) LONG RIDE > 120 km MEDIUM RIDE ≤ 120 km and > 50 km SHORT RIDE ≤ 50 km Duration Numerical (min) LONG DURATION > 300 min MEDIUM DURATION ≤ 300 min and > 150 min SHORT DURATION ≤ 150 min Average heart rate Numerical (BPM) HIGH RATE > 170 BPM MEDIUM RATE ≤ 170 BPM and > 130 BPM LOW RATE ≤ 130 BPM Intervals training? Categorical (Yes, No) INTERVAL NO INTERVAL Pains? Categorical (Yes, No) PAINS NO PAINS Motivation? Categorical (Yes, No) MOTIVATION NO MOTIVATION Sleeping? Categorical (Yes, No) GOOD SLEEPING BAD SLEEPING Eating? Categorical (Yes, No) GOOD EATING BAD EATING Table 4: Discretization of our attributes LONG RIDE MEDIUM DURATION HIGH RATE PAINS MOTIVATION GOOD SLEEPING NO INTERVAL BAD EATING Table 5: An example of one transaction more experiments we encounter that more transactions we [5] Laurel T MacKinnon. Overtraining effects on have, more rules we can obtain. In line with this, we do not immunity and performance in athletes. Immunology obtain only rules connected with overtraining, but also some and cell biology, 78(5):502–509, 2000. others tackling the basic habits of athletes, like rule 1. On [6] Shona L Halson and Asker E Jeukendrup. Does the other hand, rule 5 shows a basic example of overtraining. overtraining exist? Sports medicine, 34(14):967–981, 2004. 6. CONCLUSION [7] Lucille Lakier Smith. Overtraining, excessive exercise, and altered immunity. Sports Medicine, 33(5):347–364, In this paper we developed a prototype solution for predict- 2003. ing overtraining of an athlete. We used an Apriori data [8] Dianna Purvis, Stephen Gonsalves, and Patricia A mining method and applied it to the real dataset that was Deuster. Physiological and psychological fatigue in created by mountain biker. Experiments showed that this extreme conditions: overtraining and elite athletes. method is interesting for such manners. However there are PM&R, 2(5):442–450, 2010. also some limitations like size of the dataset. In our case we [9] Karen Birch and Keith George. Overtraining the took 50 trainings into the account to get the first picture of female athlete. Journal of Bodywork and Movement an athlete. In the future, we will also perform more tests Therapies, 3(1):24–29, 1999. with different athletes. [10] Iztok Fister, Duan Fister, and Simon Fong. Data mining in sporting activities created by sports 7. REFERENCES trackers. In Computational and Business Intelligence [1] Samo Rauter. Mass sports events as a way of life (ISCBI), 2013 International Symposium on, pages (differences between the participants in a cycling and 88–91. IEEE, 2013. a running event). Kinesiologica Slovenica, 2014. [11] Iztok Fister, Dušan Fister, Simon Fong, and Iztok [2] Kuno Hottenrott, Sebastian Ludyga, and Stephan Fister Jr. Widespread mobile devices in applications Schulze. Effects of high intensity training and for real-time drafting detection in triathlons. Journal continuous endurance training on aerobic capacity and of Emerging Technologies in Web Intelligence, body composition in recreationally active runners. 5(3):310–321, 2013. Journal of sports science & medicine, 11(3):483, 2012. [12] Iztok Fister Jr, Xin-She Yang, Iztok Fister, Janez [3] Riggs J Klika, Mark S Alderdice, John J Kvale, and Brest, and Dušan Fister. A brief review of Jay T Kearney. Efficacy of cycling training based on a nature-inspired algorithms for optimization. power field test. The Journal of Strength & Elektrotehniški Vestnik, 80(3):116–122, 2013. Conditioning Research, 21(1):265–269, 2007. [13] Simon Fong. Opportunities and Challenges of [4] Mary Black Johnson and Steven M Thiese. A review Integrating Bio-Inspired Optimization and Data of overtraining syndromeâ ˘ AŤrecognizing the signs and Mining Algorithms. In Xin-She Yang, Zhihua Cui, symptoms. Journal of athletic training, 27(4):352, Renbin Xiao, Amir Hossein Gandomi, and Mehmet 1992. Karamanoglu, editors, Swarm Intelligence and StuCoSReC Proceedings of the 2014 1 st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 51 Bio-Inspired Computation, pages 385–402. Elsevier, 2013. [14] Rakesh Agrawal, Tomasz Imieliński, and Arun Swami. Mining association rules between sets of items in large databases. In ACM SIGMOD Record, volume 22, pages 207–216. ACM, 1993. [15] Rakesh Agrawal, Ramakrishnan Srikant, et al. Fast algorithms for mining association rules. In Proc. 20th int. conf. very large data bases, VLDB, volume 1215, pages 487–499, 1994. [16] Jiawei Han, Jian Pei, and Yiwen Yin. Mining frequent patterns without candidate generation. In ACM SIGMOD Record, volume 29, pages 1–12. ACM, 2000. [17] Sergey Brin, Rajeev Motwani, Jeffrey D Ullman, and Shalom Tsur. Dynamic itemset counting and implication rules for market basket data. In ACM SIGMOD Record, volume 26, pages 255–264. ACM, 1997. [18] Edward R Omiecinski. Alternative interest measures for mining associations in databases. Knowledge and Data Engineering, IEEE Transactions on, 15(1):57–69, 2003. [19] Craig Silverstein, Sergey Brin, and Rajeev Motwani. Beyond market baskets: Generalizing association rules to dependence rules. Data mining and knowledge discovery, 2(1):39–68, 1998. StuCoSReC Proceedings of the 2014 1 st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 52 An Improved Algorithm of DPM for Two-dimensional Barcode Tao Gao Xiao-cheng Du Iztok Fister Jr. North China Electric Power North China Electric Power University of Maribor University University Faculty of Electrical Department of Automation Department of Automation Engineering and Computer Baoding, Hebei Province, Baoding, Hebei Province, Science 071003 071003 Smetanova 17, 2000 Maribor China China Slovenia gaotao81@foxmail.com iztok.fister1@um.si ABSTRACT 2. DATA MATRIX BARCODE Data Matrix is a two-dimensional matrix code which has Presently, the more commonly used international two-dimensional many advantages like as large information density, high ca-barcode includes the Data Matrix, PDF417, QR code, etc. pacity, and small size and so on. It is widely used in indus- The Data Matrix has the minimum size of all the barcodes, trial automation, logistics, transportation and other fields. and is especially suitable for small parts logos and can be A two-dimensional barcode includes printing and DPM bar- printing on to the entity directly, so that the sign is widely codes, which depending on the differences of the application used for small objects like drugs, integrated circuits, and at backgrounds of the two-dimensional barcode. DPM is the manufacturing processes from production lines [2]. abbreviation for Direct Part Mark and has some weaknesses, like low contrast, more noise jamming, complicated back- ground, uneven illumination during the data appreciation 2.1 The structure of the data matrix 2D bar- process. The main mission of this paper was to put forward a code series of image preprocessing methods which links binaryza- tion and morphological transformation based on handheld The Data Matrix of a two-dimensional barcode is shown in device, and achieve adaptive smoothing fuzzy and adaptive Figure 1. Its symbol’s structure is composed of position morphological transform through the DPM detection. The area and data areas. The position area includes an ’L’ solid experimental results show that the method can overcome boundary and an ’L’ dotted boundary as shown in Figure problems such as too large middle gap, the uneven illumi- 2, and has a dead zone with one data unit. Position area nation and noise. is the Data Matrix of border the two-dimensional barcode, and is mainly used for limiting the physical size, orienta- tion and symbol distortion of DM. The ’L’ dotted boundary Keywords is mainly used for finite unit structure of a symbol as well, Data Matrix Barcode, Binary Image Processing, Direct Part which can also help to determine physical size and distortion. Mark As shown in Figure 3, the data area contains the symbols which to be encoded, and contain coded information, like Numbers, letters and Chinese characters according to cer- 1. INTRODUCTION tain encoding rules. The DM code is the lattice and consists Two-dimensional barcode is a neotype barcode technology of two kinds of color, a black and white combination, and based on a one-dimensional barcode, and according to cri- every black or white square with the same size is a unit of teria stores data and symbolic information on special ge- data, representing a binary 1 and a binary 0 [5]. ometries using black and white, which are distributed along the horizontal and vertical areas of two-dimensional planar space. Due to some specialties regarding significant informa- tion storage capacity, good robustness and low cost [1] two- dimensional barcoding has been applied gradually within the commercial, transportation, finance, health care and other fields. Figure 1: Data matrix of two-dimensional barcode StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 53 barcode. So it is highly convenient for developers to use this barcode. 3.1 The Decoding Process of The ’Zxing’ Open Source Code In the open ’zxing’ source code, it needs to open its build camera first, and then photograph the barcode to be iden- tified. The color model of the image obtained is a RGB color model. However, the RGB color space has some dis- advantages like: No intuitive, asymmetric, and dependency on the hardware device [3]. So, the RGB color model should be converted to YUV color model, for the reason that it is Figure 2: Position area of DM important that the luminance Y is separate from the chroma U and V in the YUV color space. The formula for the RGB color model converting into YUV color model is as follows. • Y=0.299R+0.587G+0.114B • U=-0.147R-0.289G+0.436B • V=0.615R-0.515G-0.1B Figure 3: Data area of DM 2.2 DPM barcode DPM is the abbreviation of Direct Part Mark. It was orig- inally used for machinery and electronics industrial parts, and can record lots of information about parts like produc- tion, quality inspection. The factory embosses the DPM Figure 4: The basic process of zxing source code into an image by means of an etching method like the laser etching is usage can be extended to automobile manufac- After obtaining the YUV image, it is necessary to convert turing, pharmaceutical, medical, military firearms manage- the image into a grayscale image, and then convert it into ment, and so on. So the DPM barcode is a kind of important a binary image. In addition, the histogram needs to be ob- information within the Internet of technology. The DPM tained at the same time, for then to identifing the barcode. mark is the main carrier of the two-dimensional barcode The description above is the whole barcode image processing because the two-dimensional barcode has certain character- of the open source code identification of ’zxing’. The basic istics like large coding capacity, high density, high informa- process of ’zxing’ source code is shown in Figure 4. tion security. Compared with the two-dimensional barcode printed on paper, the generating method of the DPM bar 3.2 Program flowchart code varies, like ink-jet printing. In addition, it can also be The method put forward in this paper is based on ”Z-Xing” generated by certain other methods like laser etching, strik- project which is an open source two-dimensional code recog- ing a hitting machine, and electrochemical corrosion. The nition project of Google, its main purpose is to do some pre- material of the parts with DPM barcode sculpture varies treatment to the initial two-dimensional code image. It is from cast iron, aluminum, glass, hard plastic, to wood. achieved the image binarization by referring to the binariza- tion of Yang Shu, and attaching some appropriate process 3. THE SOURCE CODE to dotted data matrix bar code processing. The size of the Barcode recognition technology has been widely used, and point module will be get by the process of blob detection. the main open source code includes ’zxing’ open source code Combined with morphological processing the image can be and ’zbar’ open source code. In this section, the improve- identified, the program flow chart shown in Figure 6. Among ment is based on the ’zxing’ open source code. In the open them, it is necessary to focus on these folders mainly with zxing source code, there are certain advantages as follows: it the names ”android”, ”camera”, ”encode”, and ”result”. The can be installed within intelligent phones, and it not only has process of the initialization program is as follows: firstly, to high speed for identification but also has a short recognition load the main activity and create the object of the Capture run-time. In addition, it can identify the one-dimensional Activity Handler in this class, and then is object starts the and two-dimensional barcodes at the same time, and can camera to realize the automatic focusing. Secondly, to cre- also search online for products according to the recognized ate Decode handler threads, and then create the object of a StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 54 Decode handler, and this object then obtains original byte briefly for the Global-Histogram-Binarizer class. In the data from the camera, all the above during the first stage. improved Global-Histogram-Binarizer class, the thresh- After obtaining the data, it parses out the two-dimensional old is calculated by the kittler algorithm. Using the kit- barcode, and then parses out the character of the barcode, tler binarization algorithm, we deal with the gray images at the same time removing the characters without analysis of the DataMatrix barcodes during binarization processing, to be handle by ’Capture Activity Handlers’. This class is and then the global threshold (T) can be obtained. The kit- called the ’decode function’ of the main activity used for tler binarization processing formula is: For pages other than completing the analyse of the characters, and the image is the first page, start at the top of the page, and continue in finally displayed on the screen. In this way, it completes the double-column format. The two columns on the last page parsing of the barcode. should be as close to equal length as possible. T = e(x, y)f (x, y)/ e(x, y) (1) x y x y In Equation (1), f (x, y) is the original gray image, e(x, y) = max{|ex|, |ey|} represents the maximum gradient. In e(x, y), ex = f (x − 1, y) − f (x + 1, y) represents the horizontal gradient, ey = (x, y − 1) − f (x, y + 1) represents the vertical gradient [4]. The relevant procedure is shown in Figure 7. Figure 6: The improved binarization process Figure 5: The basic process of zxing source code 4. BINARY PROCESSING OF DATA MATRIX CODES 4.1 The binary image processing The key to binarization is to define the boundary between black and white. The image has been transformed into gray image, with every point represented by a gray value, which is needed for defining a gray value (namely threshold). If the value of the point is greater than that value, it can be defined as white (1), otherwise, it is defined as black (0). Figure 7: Kittler algorithm 4.2 The improvement of a binarization algo- 5. EXPERIMENTS rithm The process of test uses the Huawei smartphone as a carrier. As shown in Figure 6, in ’zxing’ barcode, the binarization The main parameters of the phone are as follows: the CPU is processing is implemented in Binarizer class, where the one- mediatek MT6592M, the frequency is 1433 MHZ, memory is dimensional code uses the get-BlackRow method, and the 2 GB. We deal with Data Matrix barcode on different spec- two-dimensional code uses getBlackMatrix method. There ifications, color, clarity, and 34 DPM pictures are tested. In are two classes generated from Binarizer class: Global-this section, we identified the Data Matrix barcode by soft Histogram-Binarizer and Hybrid-Binarizer. The im- wares created from the ’Z-xing’ package which have been im- plementations of the getBlackMatrix method for these two proved and haven’t been improved respectively, and compare classes are different, so we modify the binarization process the recognized effects with each other. The results are shown StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 55 Software Original Improved Number of images 34 34 Number of successful images 19 25 Success rate 55.9% 73.5% Number of images with the brightness difference 17 17 Number of successful images 7 11 Success rates 41.2% 64.7% Number of DM barcode images (2*2) 5 5 Number of successful images 2 4 Success rates 40% 80% Number of images sampled from artifacts directly 12 12 Number of successful images 7 9 Success rate 58.3% 75% Table 1: Comparison in Table 1. Because the method will be used in handheld de- 7. ACKNOWLEDGEMENTS vice finally, so the method need have more instantaneity. So This work is supported by National Natural Science Foun- in the experiment, shoot the image from different angles for dation of China (No. 71102174 and No. 51306058), the 100 times, then count average correct recognition rate and Fundamental Research Funds for the Central Universities speed. Through analyzing the experimental result, it can (No. 2014QN46). be known that the elapsed time is relevant to the contrast, roughness and illumination of image. The method sacrifices 8. REFERENCES efficiency for accuracy because of using improved Bernsen [1] C.-H. Cheung and L.-M. Po. Novel method. The elapsed time is acceptable through compar- cross-diamond-hexagonal search algorithms for fast ing the time with the elapsed time of standard DataMatrix. block motion estimation. Multimedia, IEEE Counting the time that each part spends through simula- Transactions on, 7(1):16–22, 2005. tion, it can be reached that the improved Bernsen method [2] T. Jinzhao and Q. Jie. Fast and precise detection of takes more time when the spot detection takes second place. straight line with hough transform. Journal of Image The size of image does not much affect the result because and Graphics, 13(2):234–237, 2008. the program adds size unitizing method. [3] Q. Q. Ruan and Y. Z. Ruan. Digital image processing (second edition). Electronic Industry Press, Beijing, 6. DISCUSSION AND CONCLUSION 2005. It can be seen from Table 1, that after improving ’zxing’ [4] S. Yang and Z.-h. Shang. A new binarization algorithm binarization, the success rates of various types of Data Ma- for 2d bar code image. Journal of Kunming University trix barcode recognition were improved, especially in the of Science and Technology (Science and Technology), test with the brightness difference, and the effect was very 1:011, 2008. obvious, which embodies the characteristics of the kittler al- [5] Y.-x. Zou and G.-b. Yang. Research on image gorithm. It can also more effectively find the area of uneven pre-processing for data matrix 2d barcode decoder. illumination. The improved algorithm is more stable and Computer Engineering and Applications, 34, 2009. adaptive, so as to improve the success rate of recognition. In addition, the improved ’zxing’ still has some problems: • It can be seen from Table 1, that the differences be- tween the success rates are higher, and the difference is around 20 percents, some even higher, and this is due to the pictures being a slightly less. Table 1 can reflect the effect of improvement, but it is not as great as the numbers reflected in the Table. • In the identification process of 2 X 2 DM code or 1 X n DM code, the recognition rate had improved after the improvement, but certain process of image recognition took a long time, so this kind of barcode recognition effect of the improved ’zxing’ is still unsatisfactory. • When we recognized the pictures that sampled from the artifacts directly, we needed to extract the tested sections manually, otherwise, the recognition effect was bad, in other word, and it was easily affected by the background. StuCoSReC Proceedings of the 2014 1st Student Computer Science Research Conference Ljubljana, Slovenia, 7 October 56 StuCoSReC University of Primorska Press www.hippocampus.si isbn 978-961-6963-03-9 Not for resale 9 789616 963039 Document Outline Fister jr., Iztok, and Andrej Brodnik (eds.). StuCoSReC. Proceedings of the 1st Student Computer Science Research Conference. Koper: University of Primorska Press, 2014 (Front Cover). Colophone Preface Contents David Jesenko, Domen Mongus, Borut Žalik ■ Spatially Embedded Complex Network Estimation Using Fractal Dimension László Hajdu, Miklós Krész, Attila Tóth ■ Graph Coloring Based Heuristic for Driver Rostering Viktor Árgilán, Balázs Dávid ■ A New Heuristic Approach for the Vehicle and Driver Scheduling Problems Jani Dugonik, Iztok Fister ■ Multi-population Firefly Algorithm Andrej Bukošek ■ Monte Carlo Path Tracing with OpenCL Pavel Fičur ■ A Garlic Clove Direction Detection Based on Pixel Gregor Vohl, Janez Brest ■ Simple Approach to Finding Repeated Patterns in Opponents Texas Hold’em No-limit Game Uroš Mlakar ■ Histogram of Oriented Gradients Parameter Optimization for Facial Expression Recognition Goran Hrovat, Milan Ojsteršek ■ Constructing Domain-specific Semantic Dictionaries to Supplement Domain-specific Knowledge Bases Iztok Fister jr., Goran Hrovat, Samo Rauter, Iztok Fister ■ Am I overtraining? A Novel Data Mining Approach for Avoiding Overtraining Tao Gao, Xiao-cheng Du, Iztok Fister jr. ■ An Improved Algorithm of DPM for Two-dimensional Barcode