ORGANIZING WIRELESS NETWORK AT MARSHALLING YARDS USING THE BEE METHOD

Dep. «Electronic Computing Machines», Dnipro National University of Railway Transport named after Academician V. Lazaryan, Lazaryana St., 2, Dnipro, Ukraine, 49010, tel. +38 (056) 373 15 89, e-mail viknikpakh@gmail.com, ORCID 0000-0002-0022-099X Dep. «Electronic Computing Machines», Dnipro National University of Railway Transport named after Academician V. Lazaryan, Lazaryana St., 2, Dnipro, Ukraine, 49010, tel. +38 (056) 373 15 89, e-mail olkdi@outlook.com, ORCID 0000-0002-7134-9416

Introduction
A wireless network is a connection of computers and other devices for information exchange without the use of wires; the connection is made through radio channels. The installation of a wireless network is necessary when the deployment of a cable system is impossible or economically impractical. Wireless networks have significant advantages over wired ones [5,[11][12]: they allow deploying the network in places where wires cannot be used and maintain a sufficient data transfer rate. Wireless networks also provide easy connection, user mobility, quick troubleshooting, and network equipment availability.
According to their size, the wireless networks are divided into ( Fig. 1) [5,[11][12]  According to the topology, wireless local area networks are divided into ad-hoc temporary networks (Independent Basic Service Set, IBSS); Basic Service Set (BSS); Extended Service Set (ESS). It is known that there are two main methodologies for WLAN deployment: wireless LAN with maximum service area; wireless LAN with maximum capacity [5,[11][12].
Scientists such as Ye. S. Skakov, V. M. Malysh, dealt with the planning of wireless networks [1]. To solve this problem, it is possible to use intelligent methods of multiagent optimization, which include bee, ant, particle swarm, bat, fireflies and other algorithms [1,2,18]. The work [1] determines that the bee method is the most suitable for solving the problem of locating base stations of a wireless network.
The bee method has some modifications [7][8], one of which is BCOi (Bee Colony Optimization based on the improvement concept). Its peculiarity is that the work with complete solutions of the optimization problem is considered, not with partial ones, as in the classical BCO method [1].
Since 2008, the GSM-R (Global System for Mobile Communications -Railway) standard has been used in European railways to ensure secure wireless communication between railway services and trains [3,6,7,16]. The GSM-R network is used to provide the European Train Control System (ETCS) [16]. In the UK, the GSM-R network has been fully operational since 2015 [16]. Today in Europe we are talking about transition to more modern technologies such as 4G/LTE and 5G [19]. In the future, the following new services are expected at railway transport: on-board and roadside HD video surveillance, multimedia dispatcher video streams, sensor information from railway infrastructure, including bridges, viaducts, tunnels, track defects, etc. Using infrared, sound and temperature sensors, information is collected and sent to a computation center [4].
At the present stage in Ukraine, research is being conducted on the information and telecommunication system (ITS) of railway transport using artificial intelligence methods [13][14][15]20]: neural and fuzzy networks, ant and genetic methods. But the results of the introduction of wireless networks at railway transport are not sufficiently presented yet. The authors considered the possibility of using WLAN technology for the implementation of mobile communication at railway stations, made a base station state chart in the Distributed Coordination Function (DCF), which is mandatory and based on the protocol, providing Carrier Sense Multiple Access with Collision Avoidance (CSMA/CA). Now it is expedient to conduct a study of the optimal number of base stations of the wireless network, in particular at the marshalling yard, which reaches the lower level (line enterprises) of the ITS of the railway transport of Ukraine.

Purpose
In our work, we plan to investigate the deployment of a wireless network and search for the optimal number of base stations according to the bee colony algorithm for marshalling yards of different capacities. Phyton language was chosen for software implementation.

Methodology
Problem statement. The marshalling yard (Fig. 2) includes arrival park, sorting park, hump and departure park. The most important part of the technological process is the hump, the capacity of which (small, medium, high and increased) depends on the number of sorting sets (2, 3, 4, 6-8, respectively). At the object under study there are various technological areas (point section, speed measuring and cut acceleration in the section of 1-2 BP (brake position) and 3-4 BP, measuring cut parameters, uncoupling control, cut mass determination) with various ground equipment: track sensors and photodetectors; track circuits; weight measuring devices; radar velocity meters and others. Let L be the total number of clients to be connected to the wireless base station at the marshalling yard. Known Mare the candidate places where WLAN base stations can be installed. In addition, all clients must be connected to base stations. Let us introduce the notation () k i BS rthe i-th base station of the wireless network with a coverage radius r, to which k clients are connected, where Let us consider the following function as a target one: moreover, it is necessary to comply with the restrictions: According to the IBSS topology ( Fig. 3) [5,[11][12], client stations interact directly with each other without an access point. Ad Hoc mode requires minimum of equipmenta wireless adapter. This configuration does not require the creation of any network infrastructure. In this mode, each node participates in routing by sending data to other nodes, so determining which nodes send data is done dynamically based on the network connection and the routing algorithm used. Ad-Hoc mode is used mainly to create temporary networks. According to the BSS topology ( Fig. 4) [5,[11][12], network nodes interact with each other not directly, but through an access point (AP), which can play the role of a bridge for connection to external cable network. All base stations in the network are interconnected by Distribution System (DS), which can be radio or infrared waves. The ESS topology (Fig. 5) [5,[11][12] makes it possible to combine several access points, i.e. connects several BSS networks. In this case, the access points can interact with each other. This topology is useful when you need to connect multiple users to a single network or connect multiple wired or wireless networks. It is advisable to use the ESS topology at the marshalling yard, as you can use several access points that interact with each other. In addition, there are wired networks at the marshalling yard, with which it is necessary to exchange information.
WLAN focused on the service area are developed taking into account the provision of maximum coverage with the minimum possible number of access points, Fig. 6. [5,[11][12]. In a typical service area oriented network, the ratio of the number of users to the number of access points is 25:1. Some typical features of WLAN are focused on the maximum service area: the use of bursty applications with low packet rate, for example, those generating requests to databases; the need for low spectral bandwidth requirements, so that the data flow rate can be reduced to a minimum of 1 and 2 Mbps; ensuring easy maintenance, as WLAN maintenance staff is small. In service area oriented networks, typical programs have low packet rate and low spectral bandwidth requirements. This approach allows many users to immediately access WLAN services while maintaining sufficient characteristics of the latter. Such options are common for small or medium-sized affiliates, when WLAN is chosen as an alternative to wired Ethernet. Easy-to-deploy WLAN provide basic LAN connections that are not required for file and printer sharing. Each WLAN access point serves approximately 25-30 users [5,[11][12].
Wireless LAN that focus on high bandwidth (Fig. 7) should provide maximum performance and packet rate for each BSS client. The cell size of the bandwidth-oriented WLAN is smaller than that of the WLAN, which focuses on the maximum service area, respectively, the location density of access points is higher. Highbandwidth-oriented WLAN are required when one uses applications that require high packet rates; uses applications that are sensitive to delays; deploys subnets of smaller scales (or several subnets in one service area); there is a high density of user placement. In such networks, the number of access points is several times greater than that for WLAN oriented on service area. The service area of each access point is much smaller than when building a network focused on the maximum service area. Each access point serves about 12 users [5,[11][12].
For a marshalling yard, it is advisable to use WLAN deployment methodologies that focus on maximum bandwidth, wireless LAN should provide maximum performance and packet rate for each client.
Choice of organization algorithm of wireless network at marshalling yard. Nowadays, there are many different swarm algorithms, which can be divided into [2,18]: swarm algorithms based on the behavior of insects and animals; swarm algo-rithms based on bacterial behavior and inanimate nature. The first include the following algorithms: particle swarm algorithm; ant algorithm; bee algorithm; fireflies algorithms; cuckoo search algorithm; bat algorithm. The second include gravitational search algorithm; intelligent drip algorithm; stochastic diffusion search; bacterial optimization.
Particle swarm algorithm is proposed to optimize continuous nonlinear functions. Ant algorithm is one of the most efficient polynomial algorithms for finding approximate solutions to route search problems on graphs. Bee algorithm is a relatively new algorithm for finding global extremums of complex multidimensional functions. According to the bee algorithm, bees explore areas that are near elite ones, which allows you to approximate the solution to optimal. The advantage of the algorithm is the possibility of efficient division into parallel processes and high speed [2].
Biological bases of the bee method. Bees in nature search for food by studying the space around their hive (Fig. 8) [17]. Review of the sources [1][2][8][9][10] showed that, as a rule, at the initial stage, several scout bees study the environment. After the search, scout bees return to the hive and inform the other members of the swarm about the location, quantity and quality of available food sources that they have found. The information exchange takes place with the help of dance on a specially designated area. If a bee watching a skate dance decides to leave the hive and collect nectar, it will follow one of the scouts to one of the previously found food sources. Such a bee becomes a busy forager. It collects nectar, while clarifying information about the amount of nectar near the source found. After collection, the forager returns to the hive and leaves the nectar there. It can then do one of the following: become a free forager, leaving its current source of nectar; continue extracting nectar from its source without recruiting free bees by dancing; continue extracting nectar from its source, while recruiting free bees. Described process continues anabated, while the hive accumulates nectar and explores new areas with potential food sources.
Bee colony algorithm, the enlarged scheme of which is presented in Fig. 9 [10]. First the colony is initialized, then the solution is searched, after which the bees exchange information in the hive with the help of dance, then the colony is initialized again until the stopping criterion is found, after which the optimal solution is found. Fig. 9. Block diagram of the bee colony algorithm Software implementation of the bee colony algorithm is performed in Python using standard libraries: Osis responsible for interaction with the OS; Sysis responsible for system functions; Randomis responsible for generating random numbers; Mathis responsible for mathematical operations; Datetimeis responsible for time conversion; Matplotlibis responsible for plotting.
The structure of the LocBS-BeeCol software model includes the following user classes (Fig. 10): Logused to debug the program and display the results, it is the parent for all classes; Field is a class for presenting a marshalling yard, where the search for base station location is performed; Hive is a hive class where bee-scouts exchange information and choose the optimal solution, the parent class is Field; Bee is a class of scout bees looking for coordinates of base stations, the parent class is Hive. The location of base stations (BS) at the marshaling yard is based on the algorithm, which is the basis of the LocBS-BeeCol software model and the scheme of which is shown in Fig. 11. First, all classes are initialized. This is followed by a cycle with the premise: whether all customers are covered by base stations. If not, then it is performed the clearing of extremum values, detection of search area coordinates for each bee, then the search for clients in a certain area and the coordinates of the base station are determined, at the end of the cycle the found clients are deleted from the list of free clients. If all clients are covered by base stations, the result will be displayed as a dictionary, the volume of which is the number of BS (res[bs_count]), and consists of structures whose fields are the location of the BS (new_bs_location), the number of clients (clients_in_area) covered by this BS, and numbers of these clients with coordinates (clients_in_area_list) in the form of dictionaries. Fig. 12 shows a corresponding fragment of the program. In the For cycle (line 28) the required number of bees is created, according to the initial conditions. After creating everything necessary for the work of program, the implementation of the algorithm itself begins. Line 30 starts the While cycle until all clients are connected to the base stations. Inside, the For cycle (line 32) takes place, in which the extremums and positions for each of the bees are reset to zero (lines 33-35). After specifying the search field (line 36), each bee makes a given number of attempts to find the optimal solution for the location of the base station (line 37). As a result of the cycle iteration, the global extremum, where it is more expedient to install the base station can be found, in this case it will be installed (line 40). Next, the client numbers that are connected to the BS (lines 41-42) are deleted from the general list, and the cycle returns to the beginning.   Determining the number of BS at the marshalling yard of high capacity. Input parameters: number of bees (bees_count) = 10; maximum number of attempts to find the optimal solution (max_retries_count) = 10; number of clients (cli-ents_count) = 350; maximum number of clients covered by one WLAN base station (bs_max_clients_count) = 12; coverage radius (bs_area_radius) = 50; field width (field_width) = 3 000; field height (field_height) = 600. The results of the model work are shown in Fig. 15.

Originality and practical value
Investigation of algorithm parameters. The quality of the obtained solutions largely depends on the choice of algorithm parameters. Therefore, a study was conducted to identify the dependences of the optimization results (the number of BS of wire network) and the time to find the optimal solution for different numbers of bees and attempts for marshalling yards (MS) of different capacity. The results are presented in Fig. 16-19. Fig. 16-17 show that increase in the number of bees (from 10 to 50) and the number of attempts to find the optimal solution by one bee (from 10 to 50) leads to a decrease in the number of WLAN base stations by an average of 6.5 and 9.3% respectively. Increasing the number of bees (from 10 to 50) reduces the search time for the optimal solution by bees by an average of 1.8 times (Fig. 18), while increasing the number of attempts to find the optimal solution for one bee (from 10 to 50 ) leads to an increase in the time to find the optimal solution by an average of 2.14 times (Fig. 19).
Investigation of the coverage radius of the WLAN base station was performed for a highcapacity marshalling yard. The obtained results on the LocBS-BeeCol software model are summarized in Tab. 1. The Table shows that twofold increase in the base station coverage radius (from 50 to 100 m) leads to a decrease in the number of BS approximately twice (from 136 to 64), while the search time for the solution of bees increases 2.5 times (from 8.4 to 20.6 sec). The study was conducted with the following parameters: the number of bees (bees_count) = 10; maximum number of attempts to find the optimal solution (max_retries_count) = 10; number of clients (clients_count) = 350; maximum number of clients covered by one base station (bs_max_clients_count) = 12; field width (field_width) = 3,000; field height (field_height) = 600.

Conclusions
1. Taking into account the topology of the marshalling yard and the peculiarities of the relevant technological process, it is advisable to use a wireless local area network of the ESS topology with WLAN deployment methodology, focused on maximum bandwidth.
2. Based on the results of the review of scientific sources to determine the optimal number of WLAN base stations and their location at the marshalling yard, the bee method was used, the advantages of which include the possibility of effective division into parallel processes; high speed work.
3. Based on the bee colony algorithm in Python, the corresponding LocBS-BeeCol software model was compiled, the input data of which are as follows: marshalling yard parameters (area, number of clients to be connected to base stations); WLAN parameters (base station coverage radius; number of clients served by the base station); the bee colony algorithm parameters (number of bees, number of attempts to find the optimal solution by the bee). For example, to connect 300 clients to a mediumcapacity marshalling yard with an area of 2 500x500 m 2 , 93 base stations with a coverage radius of 50 m are required.
4. The quality of the obtained solutions largely depends on the choice of the bee colony algorithm parameters. A study of the number of WLAN base stations (solution search time by the bees) on the different number of bees (attempts to find the optimal solution by bees) for marshalling yards of different capacity. It is determined that increase in the number of bees and the number of attempts to find the optimal solution by the bee leads to an increase in the quality of the optimal solutiondecrease in the number of base stations. In addition, increase in the number of bees leads to a decrease in the time to find the optimal solution and increase in the number of attempts will lead to an increase in the time to find the optimal solution by bees.
5. For a high-capacity marshalling yard in the case of increasing the coverage radius of the WLAN base station twice the number of BS is reduced by about two times, while the time to find the optimal solution for bees increases 2.5 times.


Research on Computer Aided Computation of Infectious Disease SIR Model Algorithm Based on Parameter Control

The algorithm proposed in this paper conceives the process of solving optimization problems as a process in which an infectious disease spreads among several individuals in an ecosystem. Its propagation law can be described by the SIR infectious disease model. Infectious diseases attack certain locations in several disease-causing genes of individuals. For different individuals, which disease genes and which sites are attacked are completely random; if an individual is cured, which immune genes and which sites are immunized are also completely random. Because of the problem of determining the growth rate of the number of infectious diseases, and improved SIR model is established, combined with historical data to quantitatively predict the future growth rate of the number of infected people, to formulate corresponding disease prevention and control measures.

Introduction
Infectious diseases seriously threaten human health. Diseases such as AIDS, smallpox, avian influenza, and tuberculosis have claimed the lives of thousands of people. It is increasingly important to study the transmission mechanism of infectious diseases and take relevant measures to control infectious diseases. Infectious diseases have been the focus of research in recent years, and basic theories have been applied to the fields of social science, biology, and medicine. When studying infectious diseases in groups, factors such as infection rate, mortality, and immunity cannot be ignored. Infectious Disease Dynamics It is a method for quantitative research on infectious diseases. Based on population characteristics, disease occurrence, and transmission characteristics, a model is established, and the model is theoretically analyzed to predict the development and changes of the disease, which can effectively control the spread of the disease.
Literature introduces the pulse vaccination theory in the infectious disease model in detail; literature analyzes the SIR epidemic disease model with constant population and compares the different effects of constant vaccination and pulse vaccination; literature The SIR infectious disease model of pulse vaccination was analyzed in detail; literature introduced that when the treatment resources are limited and the number of patients exceeds a certain threshold, even if the guaranteed basic regeneration number is less than 1, the disease will continue. Based on the aforementioned literature, this paper chooses the SIR model as the basic model, and assumes that the treatment resources are sufficient (to avoid backward branching), establishes a corresponding warehouse model for different strategies for controlling infectious diseases, and studies basic regeneration Number and sufficient conditions for stable diseasefree equilibrium point, compare the conditions of threshold changes under different strategies, to explore the impact of different strategies on epidemic control. Through comparison, we found that continuous vaccination is more effective than a disease in treating disease, and when vaccination When the period T <Tc, pulse vaccination is more effective for disease control than continuous vaccination [1].
At present, there are mainly four kinds of research on infectious diseases: descriptive research, analytical research, experimental research, and theoretical research. Infectious disease dynamics is an important method for theoretical quantitative research on infectious diseases. It is based on the characteristics of population growth, the occurrence of diseases, and the laws of transmission and development within the population, as well as the social factors related to them, establish a mathematical model that can reflect the dynamic characteristics of infectious diseases. Through the qualitative and quantitative analysis of the dynamic behavior of the model Analysis and numerical simulation to show the development process of the disease, reveal its epidemic law, predict its changing trend, analyze the causes and key factors of the disease epidemic, seek the optimal strategy for its prevention and control, and provide a theoretical basis for people's control decisions And quantitative basis. The thesis is aimed at the number of infectious diseases to determine the growth rate, establish an improved SIR model, and combine historical data to quantitatively predict the growth rate of the number of future infections, to formulate corresponding disease prevention and control measures.

Parameter introduction
The SIR warehouse model is to divide the population in the area into the following three categories (ie three warehouses) for a certain type of infectious disease: First, the number of susceptible people is recorded as ( ) S t , indicating that it is infected at the time but may be affected by this category The number of people infected with the disease. Second, the number of infected persons is recorded as ( ) I t , indicating the number of people who have been infected as patients at the time and is infectious. Third, the number of removed persons is recorded as ( ) R t , indicating that they have been infected from the infection at the moment the number of people who moved out of the class. The Kermack-McKendrick model is based on the following three assumptions: Hypothesis 1: Set the total population to ( ) N t , then there is The population dynamics factors such as infection, death, and migration are not considered, which means that a closed environment is considered and it is assumed that the change of disease with time is more than that of infection and death. The latter can be ignored so that the total population ( ) N t of this environment always maintains a constant       Hypothesis 3: at time t, the number of people removed from the infected person per unit time is proportional to the number of patients, and the proportionality coefficient is  , so the number of people removed per unit time is Is the number of patients removed per unit time the proportion is called the removal rate coefficient, which is also referred to as the removal rate when it is not confusing. When the removed person includes the rehabilitation person, the removal rate is also called the recovery rate coefficient or simply the recovery rate [2]. Under the above assumptions, it's the process is shown in Figure 1: (1) Further hypothesis: during the epidemic period, the population's infection rate coefficient and natural mortality coefficient are equal, and its value is  . The mortality rate coefficient due to disease is  . The incidence of the disease is bilinear βSI. The recovery rate coefficient is  The infected person is no longer infected after moving out. Figure 2 shows the block diagram of the epidemic of infectious diseases.
The first two equations do not contain R, so only the first two equations need to be discussed The basic regeneration number is 3. Infectious disease SIR model with genetic variation Infectious diseases are usually caused by viruses or bacteria. The spread of infectious diseases caused by bacteria often conforms to SIS, SIRS, and other infectious disease models, while the spread of infectious diseases caused by viruses often conforms to SI, SIR, and other infectious disease models.
There have been a lot of research results. In recent years, people have gradually introduced some prevention and control measures for disease transmission into infectious disease models for analysis and research. International research on the dynamics of infectious diseases has progressed rapidly, and a large number of mathematical models have been used to analyze various infectious diseases. Most of these mathematical models apply to the study of the general laws of various infectious diseases [3].

Related parameters
The parameters introduced in the infectious disease model with genetic variation are as follows: A represents the constant input rate to the population. Q t t  represents the probability that the genetic variant of the infectious disease will still have immunity within the period t before returning to the susceptible class after losing immunity. According to the assumption, there is Where 0 ( ) V t represents the number of individuals with genetic variation of infectious diseases that have immunity from time 0 to time f. The function 0 ( ) V t is a non-negative, monotonic, non-increasing piecewise continuous function. If the function ( ) . Therefore, the integral equation (1) Corresponding to the ordinary differential equation where  is the immunization period. This step function ( ) Q t means that the individual with genetic variation of the infectious disease at the initial moment no longer has immunity after the moment  , Therefore, whent   , the integral equation (1) Where ( ) d t u e   represents the probability that an individual with genetic variation of infectious disease at time u is still alive at time t. To ensure the continuity of ( ) V t , assume Then the integral equation (7) is equivalent to the delay differential equation In system (9), the first two equations do not explicitly contain variables ( ) V t and ( ) Q t , so the dynamic behavior of (5) can be determined by the following delay differential system According to the actual background, assume that the initial condition of (9) is 1 2 It is easy to know that there is a unique solution to 0 t  under the initial conditions (8)

Stability analysis of an epidemic model with genetic variation
3.2.1. Analysis of differential systems with delay. The Yizhi system (9) has a positively invariant set . At the same time, the following conclusions can be obtained: The system (9) always has a disease-free balance point There is only one endemic disease balance point when 1 1 P is globally asymptotically stable on D when 1 1 R  ; D is globally asymptotically stable on 2 P when 1 1 R  . For the system (9), the existence of the equilibrium point can be known directly. The global stability of the disease-free equilibrium point 1 P is discussed below. 1 x S S   (15) Then the system (9) becomes ' 1 3) Define the Lyapunov functional Therefore, when 1 1 R  is the endemic disease equilibrium point 2 P is globally asymptotically stable.

Case analysis
This part takes a sudden outbreak of an infectious disease in a boarding school as an example. The empirical method of determining the parameters of the SIR model; in a boarding school in the UK, the whole school The number N=763. In an outbreak of infectious disease, more than 500 people have been infected in the school within 14 days (no deaths). Infected people will be forced to suspend classes. Because it is a boarding school, the total population is relatively stable. Meet the basic assumptions of the SIR model [4]. The implementation of the algorithm here is also based on the previous parameter assumptions in the article, the parameter vector     2 = , 0,10     , which is accurate to three decimal places when the parameter accuracy is determined; the combination of path update and parameter update adopted in the  Figure 5) are the parameter inversion results obtained when the initial population number is N=20,10. In this example, although the data dimension is low and the amount of data information is low, the GA algorithm is still A good parameter estimation result is achieved; when the initial population number is selected as N = 20, in the actual inversion process, the algorithm can obtain a stable convergence solution    , It is consistent with the basic regeneration rate of 3.7143 obtained in, and the relative error is about 0.12%. The reference value in the specific control and treatment of the disease is also basically consistent. It is considered that the GA algorithm has achieved very good parameters to Estimate the effect. The table presents the parameter estimation results obtained by the algorithm used in [6] and the algorithm used in this paper: 3.7187 As can be seen from the table, the estimation results obtained by using the GA algorithm in this paper ((1.666, 0.448) are the same (this is related to the selection of the fitness function in the text), indicating that the GA algorithm also implements the parameter estimation of the SIR model. It has certain stability and practicability; and the algorithm is convenient to get started, even in high-dimensional parameter estimation, it has considerable advantages, so this paper believes that the GA algorithm can be widely applied to similar infectious disease dynamic models Parameter estimation. In the end, it is worth affirming that any model or algorithm is just a simplified simulation of a real problem, limited to various assumptions in the modeling process. We must start from a variety of methods for comprehensive analysis of any problem, the results obtained Sexual parameters give in-depth explanations and multifaceted parameter references, and should not be rushed to be used directly for practical problem guidance [7].

Conclusion
This article has also studied a SIR infectious disease model with genetic variation. Most infectious diseases such as smallpox, flu, hepatitis, measles, etc. have strong immunity after being cured, so the people who are cured are not healthy (susceptible Person), not a patient (infected person), they have withdrawn from the Department of Infectious Diseases, and, for many epidemics, susceptible people are vaccinated, but the classic SIR model does not consider this situation. In this regard, a SIR infectious disease model with vaccination was built based on the classic SIR model, and its asymptotic stability was proved through theoretical analysis. The disease-free balance point and endemic disease balance point of the model were discussed. Through the linearization method and Lyapunov function, the following conclusions are obtained: There is a disease-free equilibrium point 1


extfieldsofstudy:
Materials Science
Medicine
Physics
Mathematics

text:
Modelling the dielectric constants of crystals using machine learning

The relative permittivity of a crystal is a fundamental property that links microscopic chemical bonding to macroscopic electromagnetic response. Multiple models, including analytical, numerical and statistical descriptions, have been made to understand and predict dielectric behaviour. Analytical models are often limited to a particular type of compounds, whereas machine learning (ML) models often lack interpretability. Here, we combine supervised ML, density functional perturbation theory, and analysis based on game theory to predict and explain the physical trends in optical dielectric constants of crystals. Two ML models, support vector regression and deep neural networks, were trained on a dataset of 1,364 dielectric constants. Shapley additive explanations (SHAP) analysis of the ML models reveals that they recover correlations described by textbook Clausius-Mossotti and Penn models, which gives confidence in their ability to describe physical behavior, while providing superior predictive power.

I. INTRODUCTION
The dielectric response function is one of the fundamental properties of material, which can give an insight into optical and electric properties. Specifically, the electronic highfrequency component of the dielectric constant has been explored intensively in early studies and is still of interest among researchers today. 1,2 Numerous efforts have been made to model dielectric constants. Most notably, the Clausius-Mossotti and Penn models, and their variants, are widely adapted throughout the literature. [2][3][4]

A. Clausius-Mossotti Model
The Clausius-Mossotti (CM) equation expresses the dielectric constant ε as: Here, α i is polarizability of atomic species i, and v is the volume. In the case of molecular crystals, α i can be assigned to a constituent molecule, [5][6][7] and in ionic solids it is assigned to an ion. The inherently many-body nature of the dielectric constant is reduced to this simple relation through employment of two large assumptions.
a) Electronic mail: keith.butler@stfc.ac.uk b) Electronic mail: a.walsh@imperial.ac.uk Firstly, that an external electric field is screened by a dielectric medium before reaching an atom. In the CM model, cubic symmetry is assumed and this local field E loc (r) is expressed as where E ext (r) is the external electric field. This relation, the Lorentz relation, is frequently used for non-cubic cases and generally holds when the anisotropy is small. 8,9 Secondly, the atomic polarizability α i is assumed to be a constant additive quantity. In other words, α i is assumed to be unaffected by the surroundings. Due to its simple form, the CM model is still being used and in practice works for many materials, 2,10 despite the underlying assumptions. 8,[11][12][13] Numerous efforts have been made to improve the CM model. [14][15][16] In general, more effort was put into improving atomic polarizability α i rather then changing the functional form. 14-16 While Wilson and Curtis increased the accuracy of the model by considering electrostatic environment of an atom, 15 Jemmer et al. fitted parameters to Møller-Plessett calculations. 16 Nonetheless, the range of compounds and structures in which the model is valid was found to be limited. More recently, Shannon and Fischer have fitted atomic polarizabilities to sets of experimental and computational data, consisting of 650 values from 487 compounds. 2 This was one of the first studies to consider compounds throughout a wider variety of structures, and their values of atomic polarizability were more transferable. However, due to the limited dataset, some parameters were manually fitted using chemical intuition. Despite the success of the CM based models for many compounds, 2,15,16 they are far from being universal and poorly describe covalent systems.

B. Penn Model
The premise of the Penn model is the electronic band picture. When the overlap between electrons at different atomic sites is large, such a model is anticipated to be favorable. Treating the external electric field as a time-dependent perturbation, Penn derived the following relation: 1 where ω p is the plasma frequency and E g is the width of the band gap. The significance lies in the fact that the dielectric constant is described only by the plasma frequency and the band gap. Using the second order perturbation theory, polarizability α for an isolated atom can be written α ≈ ( ω p /E g ) 2 , where E g is the HOMO-LUMO gap. In this case, the screening of the external field can be ignored and Penn model is equivalent to the CM model.
As with the CM model, there are various assumptions for the Penn model. The first main assumption is to consider the electrons as a free electron gas. This allows the form of the equation to be greatly simplified, but is not valid when electrons are localized. 17,18 The second assumption is to consider only the bandgap and not the density of states around the valence and conduction band edges. This flat band approximation results in a large error when the valence or conduction band has a large dispersion. 4 Various efforts have been made to improve this model. Phillips added a correction for ionic bonds and showed that dielectric constants in zinc-blende and wurtzite were reproduced. 17,18 Furthermore, Van Vechten added a scaling factor taking into account that not all valence electrons contribute to the dielectric response. 19,20 These modifications did contribute to increasing the predictive power; however, a universal model was not achieved.

C. Data driven models
The CM model approaches the problem from a molecular picture, whereas the Penn model is based on the electronic band picture of a crystal. As such, the former has difficulty with delocalized electrons and the latter has difficulty with localized electrons. In both cases, the number of model parameters were increased in order to improve transferability. In all cases, this was done empirically and the physical justification was often limited.
Rapid progress in computer technology and first-principles modelling techniques has enabled the calculation of dielectric constants for large number of materials. [21][22][23] The modern calculation of dielectric constants are typically done using density functional perturbation theory (DFPT). 3,24-28 Therefore, many studies have involved analysis of datasets of dielectric constants aiming to understand trends or predict values for new materials. For example, Han and coworkers have performed high-throughput calculation of binary and ternary inorganic compounds, and compared the relation between the band gap and the dielectric constant. 3,24 Some studies employ machine learning (ML) methods where a statistical model is trained. 25,26 For example, Umeda et al. trained different ML models on a dataset of 3,382 compounds. 25 They obtain a good agreement with DFPT calculations; however, the reasons for this agreement are not discussed.
Overall, simple analytical models are limited by their generality, whereas ML models are difficult to interpret. By studying the difference between ML and analytical models, we may recover new underlying physics and deepen our understanding of the phenomenological behavior of the dielectric response.
In this study, we use two types of ML models and investigate the reasoning behind their predictions. Firstly, a deep neural network and support vector regression model are trained on the same dataset. Next, their results are compared with DFPT calculations. Lastly, we perform Game Theoretic analysis (Shapley additive explanations) to elucidate the characteristics of the trained ML models and to investigate the reason behind their predictions. We compare the result closely with the CM and Penn models and show that ML can indeed learn underlying physical trends, while having superior predictive power. Additionally the Shapley analysis allows us to identify reasons for poor model performance in particular cases, which is important when considering how relaible predictions are.

A. Machine learning
Supervised ML models were trained on a dataset of 1,364 dielectric constants. The dataset was prepared by combining two pre-existing datasets, 27,28 and averaging the data for the overlapping materials. Materials with a band gap less than 0.5 eV were removed because small gap materials require very dense sampling of the Brillouin zone, which is difficult to realize in high-throughput calculations. 33 The dielectric constant was calculated by taking a diagonal average of the electronic part of the dielectric matrix. Since this dataset only contained minimum features, it was augmented and processed using Materials Project, 29,30 pymatgen, 32 SMACT, 31 and scikit-learn. 34 The dataset was split into 8:2 the training:test sets. The same data was used when training the two different ML models.
The first ML model was support vector regression (SVR). Before training, we analyzed the available features and removed those unimportant or similar. This step was necessary, because having too many features increases the dimensions of the manifold that the model must learn in and thus lowers the performance. 35 Specifically, we calculated feature importance using the random forest ( Fig. S1) and removed low importance features: Space group (one hot encoded), atomic species (one hot encoded), number of elements, and number of sites inside the unit cell. Next, with the criterion of r 2 > 0.90, we removed maximum oxidation state and variation in Maledung energy, which had high correlation with variation in oxidation state and minimum in Madelung potential, respectively (Fig. S2). The refined features are presented in Table I.
Four machine learning models, the random forest, gradient boost regression, kernel ridge regression, and SVR, were trained and compared ( Fig. S4 and Table S1). The former two models are based on decision trees and latter two models are based on kernel methods. Since the random forest and gradient boost regression are ensemble methods, we surmise their performance to be better, however, they were overfitting and the kernel methods were able to generalise better. The kernel ridge regression had slightly lower performance compared to SVR, therefore SVR was employed for this study. The second ML model was a deep neural network (DNN). Specifically we have adopted the DNN architecture MEGNet developed by Chen, et al. 36 MEGNet can be trained only by using the crystal structure of materials. The exact structure of the network is presented in Fig. S3. It is not trivial to uniquely express a crystalline system in the form of a vector, however MEGNet overcomes this difficultly by representing bonding networks as graphs and using the set2set algorithm to consistently treat different sized graphs. 36 Information of atomic number and bond length are encoded into the graph representation. The model was trained for 400 epochs and the layer weights of the epoch epoch with the smallest Huber loss were employed.
There are two advantages and a disadvantage of using MEGNet over SVR. The first advantage is to avoid so called "feature engineering". Feature engineering is a process of feature selection, which was necessary for SVR. Although this can be a way of inputting our domain knowledge, 21 this relies heavily on experience and intuition and may obstruct systematic improvement of an ML model. The second advantage is that it is easy to perform transfer learning. As larger datasets are available for other materials properties, it is tempting to use that information in training. Since the upper layers in DNNs are known to learn general trends, it is possible to import the layer weights of the upper layers in advance and improve the accuracy. 36 Taking into account the Penn model, and the results of feature engineering for SVR (Figs. S2 and S1), we performed transfer learning from MEGNet model Chen et al. trained on band gaps. 36 The disadvantage of DNNs is that interpretability. Specially, because we used graph representation as the input, the relation between the features and dielectric constant were difficult to extract.

B. Density functional theory calculation
In order to validate our results, we have performed DFPT calculations on a total of 24 structures. These structures were the materials with largest error in the SVR and DNN models. The calculations were performed with projector-augmented wave scheme as implemented in VASP. [37][38][39] The structures were taken from Materials Project and were calculated using the PBE functional. 29,30,40 The reciprocal space were sampled so that the spacing between the k-point was about 2π × 0.03Å −1 . The energy cut off was set to at least 600 eV and the wavefunctions were optimized to a tolerance of 10 −6 eV.

C. Shapley additive explanation
For model interrogation, we perform Shapley additive explanations (SHAP) analysis. 41 The Shapley regression value is defined as: Here, S is a subset of the features (F ), i is a particular feature of interest, f is the ML model, and x S represent values of the input features in the set S. 41 φ i describes how much the model output changes when feature i is added to the model. Therefore, it can be used to quantify feature importance. As SVR cannot take features of different size, we used the median of the feature instead of removing them. Furthermore, if we fix data x S and model f , we can show which features are responsible for the prediction given subset x S . Using the additivity approximation suggested by Lundberg and Lee, 41 we calculated SHAP values of all the features and data in the dataset.

A. ML prediction
The dielectric constants predicted by SVR are plotted against the dataset values in Fig.  1(a). In general, materials with dielectric constants above 10 had larger error. DNN performed similarly to SVR as shown in Fig. 1(b). The main difference was suppression of error in the high dielectric constant range and slight increase of error in the low dielectric constant range in the predictions from the DNN.
Performance metrics of SVR and DNN are summarized in Table II. This performance is similar to the performance of ML models in previous studies. [42][43][44] Ideally, if there is no overfitting and the ML model is trained sufficiently, metrics for training and test data should be identical, however Table II imply tendency for both models to exhibit overfitting. For SVR, different sets of hyperparameters and loss functions were trialled and the result presented here represents the best achievable performance (Table S1).We built a series of SVR models on different sized subsets of the data (Fig. S6) and find that the difference in the loss function metric decreased as we increased the dataset size, but did not converge. We therefore attribute overfitting to the small size of our dataset. Compared to SVR, the difference in metric between the test and training data were smaller for DNN, suggesting that the latter model was able to generalize better, and exhibits overfitting to a lesser extent. For test data, the DNN demonstrated higher performance for all of the metrics than the SVR (Table II).
The 10 materials with largest prediction error are listed in Tables III and IV. Out of the 10 materials, seven of them were common amongst the two ML methods. It is worth noting that LiBC, Ga 2 Te 5 , LiAsS 2 and In 2 HgTe 4 have error large enough to be identifiable in Fig. 1. We also present the energy above hull (obtained from Materials Project) to show that these materials are not artefacts from computational screening and could possibly be stable. 29,30 Since the dielectric constant dataset we use for training was derived from a highthroughput workflow, 27,28 the precision is optimized to be convergent for most materials, 45 which results in some results being poorly converged. 33 Therefore, for comparison, we have performed higher precision calculations using DFPT on these compounds and the results are presented in Table III and IV.  TABLE III. List of 10 material with largest error in support vector regression (SVR) prediction. εSVR, εDataset and εDFPT(this work) is dielectric constant from SVR prediction, training dataset and our density functional perturbation theory (DFPT) calculation, respectively. ∆ε dataset is the difference between εSVR and εDataset. ∆εDFPT is the difference between εSVR and εDFPT. ∆E hull (meV) is energy above hull per atom. The larger error found when the large dielectric constant is greater is not surprising (Fig. 1). In addition to the numerical instability in DFPT calculations, 33 there are fewer materials with large dielectric constants in the dataset, therefore the ML models were not able to fully learn the trends in the large dielectric constant range.
Comparing SVR and DNN, the performance was similar, as shown in Fig. 1 and Table  II. This is surprising since the features used for training each model were different. It suggests that both models are sufficiently capturing physical trends in the training data. If the model is generalized well, we expect them to be able to detect anomalous data in the DFPT dataset. For eight (seven) materials out of 10 materials, SVR (DNN) was actually predicts values closer to our calculations, which were performed under higher precision. This is clearly the case for Ga 2 Te 5 and LiAsS 2 . For example in the case of Ga 2 Te 5 in Table  III, the SVR prediction was 10.48, whereas the dataset value was 17.19 and our calculation was 13.75. Since small band gap materials require fine Brillouin zone sampling, 33 the standardized sampling in the high-throughput calculation setup may be insufficient. To confirm this, we calculated the dielectric constant of LiAsS 2 with different Brillouin zone sampling density. Our converged value of dielectric constant was 8.57, while 11.46 was obtained using coarse sampling. 11.46 is closer to the dataset value of 12.39 (Tables III and S3). Although this is not a direct evidence, we speculate that materials which exhibited a large difference between our calculations and the dataset were especially sensitive to the Brillouin zone sampling and as a result, made the reported dielectric constants of these materials anomalous.
B. SHAP analysis SHAP analysis applies a game theoretic approach to calculate the importance of individual input features to a given model prediction. A positive (negative) SHAP value indicates that a given feature contributes to an increase (decrease) in the prediction with respect to the mean of the set. Fig. 2 shows calculated SHAP values for all features and data. The features are ordered by their importance. Note that the order of importance differs from Fig. S1, which was a random forest model trained on all features and data, whereas Fig. 2 is SVR built only on training data.
The high contribution of material density can be explained by both the CM and Penn expressions. If there are more electrons in a given volume, the dielectric response will become larger, and indeed SHAP analysis shows that dielectric constant monotonically increases with density.
The high contribution of the band gap is also no surprise from the Penn model (Eq. 3). Lower energy excitations result a larger dielectric constant. A large band gap gives a negative SHAP contribution. Interestingly, the magnitude of positive contributions from small band gaps has a longer tail in the distribution than the negative contribution from large gaps. This suggests that although a large gap decreases the dielectric response, this effect diminishes quickly so that continuously increasing band gap will not always decrease the dielectric constant.
To take a deeper look at these relationships, further analysis was performed in Fig. 3(a). The contribution from band gap decreases as the gap increases with an inverse power law relation. By coloring the individual points according to material density, the interplay of density and band gap can be observed. Specifically when the band gap is low, a low density increases the SHAP value, while when the gap is high a low density decreases the magnitude of the negative SHAP value. While it is not possible to derive rigorous analytical relationship from number of data points we have, Fig. 3(a) suggests a φ Eg ∼ n/E m g relation, where φ Eg is the SHAP value for the band gap, n is the material density, and m is an arbitrary constant. This is similar to the Penn relation (Eq. 3). Therefore, it is possible to interpret that ML is learning the Penn model while incorporating other feature relations as correction terms.
When the formation energy is high (low) the dielectric constant is raised (suppressed),  Table III. Starting with the dielectric constant predicted by a median of the features (ε = 5.58), each feature contributes to lowering or elevating the value of the dielectric constant in an additive manner. The final value of the dielectric constant ε pred is the output of the model. Fig. 2. This relation is opposite to that of the band gap and the negative correlation −0.65 between the two also agree with this trend (Fig S2). Within our ML model, two roles of formation energy can be suggested. Firstly, it is acting as a ionicity parameter. According to Pauling, the formation energy of binary compounds is ∆H ∝ (∆E Pauling ) 2 , where ∆E Pauling is the Pauling energy difference. 46 This trend was present in our dataset (Fig. S7); however, the variance was large which suggests contributions from other features. Secondly, that it is another feature for the band gap. In main group binary compounds, the relation ∆H ∝ (E 2 g /E f ) ln(E g /E f ) has been reported. 47 Quantitatively, this relation did not hold for our case (Fig. S7). Given these two relations and accounting for the fact that formation energy did not show obvious relation between the Pauling energy difference and the band gap, we suggest that it contains both information weakly.
The variation of oxidation states also exhibited high importance, Fig. 2. Since oxidation states have information about ionicity, we may expect it to lower the value of the dielectric constant when it is large. 19,20 Fig. 2 shows the opposite trend. This counter intuitive result can be explained by interpreting this feature as a correction the flat band approximation, described earlier for the Penn model. As the band gap is defined as a difference between the valence band maximum and the conduction band minimum, it does not have information about band dispersion. If ionic and covalent compounds with the same band gap exist, ionic compounds will have smaller effective band gap. When band gap is large (small) and contributes to making dielectric constant smaller (larger), lower variation in oxidation states amplifies (suppresses) this effect (Fig. S8). Theoretical studies also support the importance of this correction. 4 Other oxidation state and Madelung energy features did not show clear trends, which constitutes further evidence that the variation of oxidation state is acting as an ionicity parameter. This result would not have been available if we were only considering feature importance, and highlights the ability of SHAP to access deeper understanding of trained models.
Since the distribution of ionic species ratio features had broad SHAP value distributions and were difficult to interpret in Fig. 2, we continue their discussion based on Figs. 3(bd). SHAP values of O 2and Fratio are plotted against material density. In general, higher concentration of these ionic species tend to reduce the dielectric constant. This effect could be interpreted as originating from the strongly electronegative nature of O 2and F -, which thus stabilize valence electrons and suppress their dielectric response. In contrast, for the Li + cation, the SHAP value is positive, meaning it increases the dielectric constant. Furthermore, we can see the interplay of density and composition; when the density is higher, the effect of O 2and Fratio was enhanced, where O 2had slightly larger change. Due to the limited number of data points, further study is required to draw firmer conclusions.
Finally, we explore the limitations of the trained SVR model by examining LiBC and CdCN 2 in Table III. Since SHAP values are additive, 41 we are able to decompose them to each of the features. We begin with the dielectric constant predicted for the median of features (ε = 5.58) and add contribution of all the features sequentially (Fig. 4). Although only the 10 most important features are presented, other features had minor contribution ∼ 0.001, so we will only discuss features presented in Fig. III. For LiBC, which had largest error for both SVR and DNN, SVR predicted 3.54 where our calculated value was 11.96 (Table III). As shown in Fig. 4(a) the largest contribution was from the material density. Together with the fact that the material density was small (2.15 g/cm 2 ), the error is likely due to the incorrect contribution of material density. It is likely that in this range of material density, sufficient sampling was not achieved with our training dataset.
Most of the materials with large error feature a large contribution from material density. The exception was CdCN 2 . SVR predicted 6.55, where 4.45 was the result of our DFPT calculation (Table III). As shown in Fig. 4(b), the three largest contributions are N 3ratio, formation energy, and variation in oxidation state, which have reasonable values. We instead assign the error as arising from the large anisotropy, which is not taken into account. The full dielectric matrix is: The DNN performed better because the bond lengths are explicitly taken into account through the graph representation.

V. CONCLUSIONS
We showed that two different machine learning models, a support vector machine and a deep neural network, were able to predict the dielectric constants of crystals with reasonable precision. Comparison with our density functional perturbation theory calculations reveals that the machine learning models were able to detect erroneous results in the original dataset. We performed SHAP analysis of the support vector machine model, which illustrated that it is learning similar relations to the textbook Clausius-Mossotti and Penn expressions. Finally, we showed the limitation of our support vector model through detailed analysis of the predictions.
We suggest that as long as the dataset is sufficiently large to sample the crystal space of interest, machine learning models can be an effective approach not only for prediction of material properties but also to capture physical trends. The analysis approach used in this study is not restricted to dielectric response, and therefore has application potential for other relations, including properties that are not as intensively studied and where there is no existing analytical description.

VI. DATA AVAILABILITY STATEMENT
A set of electronic notebooks to reproduce the model training and analysis performed in this study are available from [DOI to be inserted].


Contextual-Relation Consistent Domain Adaptation for Semantic Segmentation

Recent advances in unsupervised domain adaptation for semantic segmentation have shown great potentials to relieve the demand of expensive per-pixel annotations. However, most existing works address the domain discrepancy by aligning the data distributions of two domains at a global image level whereas the local consistencies are largely neglected. This paper presents an innovative local contextual-relation consistent domain adaptation (CrCDA) technique that aims to achieve local-level consistencies during the global-level alignment. The idea is to take a closer look at region-wise feature representations and align them for local-level consistencies. Specifically, CrCDA learns and enforces the prototypical local contextual-relations explicitly in the feature space of a labelled source domain while transferring them to an unlabelled target domain via backpropagation-based adversarial learning. An adaptive entropy max-min adversarial learning scheme is designed to optimally align these hundreds of local contextual-relations across domain without requiring discriminator or extra computation overhead. The proposed CrCDA has been evaluated extensively over two challenging domain adaptive segmentation tasks (e.g., GTA5 to Cityscapes and SYNTHIA to Cityscapes), and experiments demonstrate its superior segmentation performance as compared with state-of-the-art methods.

Introduction
Semantic segmentation has been a longstanding challenge in computer vision, which aims to assign class labels to every pixel of an image [59]. Deep learning based approaches have achieved great successes at the price of large-scale densely-annotated datasets [3,27,9] which are prohibitively expensive to collect [9]. One way of circumventing this constraint is to use synthesized images with automatically generated labels (e.g., synthesized [36] or game-engine produced  [44,49,45] (e.g., probability-/entropy-/patch-represented output), which tends to introduce segmentation errors due to the neglect of local contextual consistency. Our CrCDA shown in the lower part adapts features at local level for contextual-relation consistency between the source and target domains which generates more accurate segmentation consistently. In the graph, "compatibility relations?refer to visual patterns with high co-occurrence frequency (e.g., "pole?should be beside the "sidewalk"), and "incompatibility relations?refer to visual patterns with low cooccurrence frequency (i.e., "sky" should not in the "building").
[35] data) in network training. Unfortunately, such models usually undergo a drastic performance drop when applied to real-world images [53] due to the domain bias and shift [39,32,44,48,40,30]. Unsupervised domain adaptation (UDA) has been introduced to address the domain bias/shift issue. To reduce the cross-domain discrepancy, most state-ofthe-art UDA methods [46,44,31,49,45,17] exploit adversarial learning for distri-bution alignment in the intermediate feature [46], output [44,31] or latent [45,49] space. Among this cohort of adversarial-based methods, a common and pivotal step is the employment of a discriminator [16] that predicts a domain label for data being either source or target domain. However, the discriminator works only on image-level and merely achieves global consistency (i.e., locational/spatial distributions consistency), where local contextual consistency (i.e., region-wise contextual-relationships) is largely neglected.
Local contextual-relationships are ubiquitous and provide important cues for scene segmentation. They can be formulated in terms of semantic compatibility/incompatibility relations between one thing/stuff and its neighbouring things/stuff. Under this formulation, a compatibility relation is an indication of visual patterns with high co-occurrence frequency, e.g. a pole beside a sidewalk, and an incompatibility relation is an indication of visual patterns with low co-occurrence frequency, e.g. a person above a driving car. The contextual information has been extensively explored in supervised semantic segmentation, whereas the local contextual-relationships is largely neglected in unsupervised domain adaptive semantic segmentation though theyre beneficial in addressing local contextual consistency and inconsistency in the target domain, as illustrated in Fig. 1.
To this end, we propose an unsupervised domain adaptation method for semantic segmentation that explicitly models the local contextual-relations in the feature space of source domain (with label) and then transfers this contextual information into the target domain (without label), ultimately improving target domain segmentation quality, as shown in Fig. 1. We first establish local contextual-relationships pseudo annotations in the source domain. This can be achieved by sampling regions from pixel-level ground-truth maps of source images and clustering the sampled regions to indexed N /M groups via Dbscan [12], as illustrated in Fig. 4. With the local contextual-relationships pseudo annotations in source domain, we can train a classifier C to explicitly models/learns the local contextual-relations in the feature space of source domain, and then transfers/enforces these local contextual-relations into target domain.
Following current discriminator-based global alignment methods [44,45,49,31], a intuitive idea is to employ hundreds of discriminators to align hundreds of contextual-relations across domain where a single discriminator focuses on a single contextual-relation, or employ just one discriminator to align all contextualrelations across domains. Obviously, the former is cumbersome which requires much redundant computation, while the latter is not aware of a variety of contextual-relations in the data distribution and may end up biasing to lowlevel/simple difference. Therefore, different from current discriminator-based global alignment methods [44,45,49,31], we enforce these local contextual-relations on target domain via adaptive entropy max-minimizing (AEMM) between classifier C and feature extractor E that estimates prototypical feature representations of these local contextual-relations and congregates neighboring target incorrect samples/contextual-relations to the approximated correct source prototypes alternatively, ultimately leading to consistent local contextual-relations across do-mains. In this way, our method requires no discriminator which is normally used in UDA-based semantic segmentation and introduces training instability and extra components. In addition, this AEMM learning scheme can also be applied into pixel-/global-scale training.
The contributions of this work can be summarized in three aspects. First, we propose an unsupervised domain adaptation method for semantic segmentation that explicitly models the local contextual-relations in the feature space of source domain (with label) and then transfers this contextual information into target domain (without label). To the best of our knowledge, this is the first effort to explore contextual information for UDA-based semantic segmentation. Second, it introduces a novel adaptive entropy max-minimizing adversarial learning scheme to effectively align hundreds of local contextual-relations across domain, which requires no discriminator and adds no overhead. Third, it shows the proposed method can be seamlessly integrated into existing domain adaptation techniques without extra overhead except two classifiers and achieves consistent improvements on semantic segmentation. Fourth, extensive evaluations over two challenging UDA tasks GTA5 → Cityscapes and SYNTHIA → Cityscapes show that our method achieves superior semantic segmentation performance consistently.
Adversarial learning based approach: Adversarial learning based UDA has been extensively explored for semantic segmentation, where a discriminator is employed to minimize the divergences between source and target domains in feature or output spaces. [19] first applies adversarial learning for UDA based semantic segmentation by aligning feature space at global scale. Curriculum domain adaptation [55] utilizes certain inferred properties (e.g., superpixel and global label distributions) as the guidance to train the segmentation network. In [44] and [7], the adversarial learning is used to align the global structure to benefit from the scene layout consistency across domains, where [7] integrates a target guided distillation module to achieve style adaptation. In addition, [38,39] combines adversarial learning and co-training to achieve domain adaptation via maximizing the discrepancy between two classifiers' outputs.
Image translation based approach: Inspired by the recent advances in image synthesis (e.g., CycleGAN [60]), a number of GAN-based methods are proposed to generate target images conditioned on the source, which can help reduce the domain discrepancy before training segmentation models. CyCADA [18] uses CycleGAN to generate target images conditioned on the source images and achieves input space adaptation with a joint adversarial learning for feature alignment. A similar method, DCAN [50], implements channel-wise feature alignment to preserve spatial structures and semantic concepts in the generator and segmentation network. [42] transfers the information of the target domain to the learned embedding via the joint adversarial learning between generator and discriminator. Besides using GANs [16] to align the embedding across domains, [62] proposes a novel conservative loss to penalize the extremely easy and difficult cases while enhancing moderate examples.
Re-training based approach: Another approach of UDA based semantic segmentation is pseudo label re-training [64,63,26] that uses high-confident predictions as pseudo ground truth for the target unlabelled data to finetune the model trained on the source data. In [64], class balancing and spatial prior are included to guide the iterative re-training in target domain. [49] proposes a soft-assigned version of re-training, where it enforces the "most-confused" pixels (e.g., with equal probabilities for all classes) to become more confident (i.e., with either low or high probability for each class) by entropy minimization. [64] instead implements iterative learning on high-confident pixels.
Our method does not follow either global/class-wise feature space alignment using discriminators [19,28,29,46,31,7] or re-training on target data [64,41]. Instead, we enforce multi-scale feature space alignment via multi-scale entropy max-minimizing. To the best of our knowledge, this is the first end-to-end multiscale UDA network that achieves competitive performance on two challenging UDA tasks.

Methods
In this section, we present our framework for contextual-relationships consistent domain adaptation (CrCDA): a discriminator-free adversarial training scheme between a feature extractor module and a classifier via adaptive entropy maxminimizing (AEMM) to align local contextual-relationships across domains. Fig.  2 illustrates our network architecture.

Problem Definition
We focus on the problem of unsupervised domain adaptation (UDA) in semantic segmentation. Given the source data X s ⊂ R H×W ×3 with C-class pixel-scale segmentation labels Y s ⊂ (1, C) H×W (e.g., stimulated images from game engines) and the target data X t ⊂ R H×W ×3 without labels (i.e., real images), our goal is to learn a semantic segmentation model G that performs well on the target dataset X t . Current adversarial learning methods rely heavily on discriminators to align the distributions of source and target domains via two loss functions: segmentation loss on source data and adversarial loss for alignment.
However, there exists a crucial limitation for these approaches: even if perfect adaptation is achieved through a discriminator, the alignment is implemented on global level (i.e., image-level), where local contextual information may be lost/deconstructed. The reason lies in that the discriminator can only implement alignment at global level, which inputs the whole map but outputs a digit to represent domain labels (e.g., 0 or 1). In some cases, parts of local regions (i.e., local contextual-relations) have been well aligned across domains. However, the discriminator might deconstruct this existing local alignment during implementing the global marginal distribution alignment. In this paper, we define this phenomenon as "lack of local consistency (i.e., local contextual inconsistency)", which is important to semantic segmentation in dense pixel-scale prediction.

Overview of Network Architecture
As shown in Fig. 2, our semantic segmentation model G consists of a feature extractor E and two classifiers (i.e., C seg and C cr ) where C seg is for pixel-scale segmentation and C cr is for local-scale contextual-relations learning/classification. E extracts features from input images. C seg and C cr classify features generated by E into pre-defined semantic classes. Specifically, C seg processes features at pixel-scale, which aims to predict pixel-scale labels. The pre-defined semantic class domain for C seg is the pixel-scale ground-truth, so there is no difference between C seg and traditional segmentation classifier. C cr processes features at local scales, which aims to predict region-scale/contextual-relations labels. The pre-defined semantic class domain for C cr is the clustered contextual-relations ground-truth. The establishment procedure of clustered contextual-relations labels is described in Section 3.3 and shown in Fig. 3. We train E and the classifiers (i.e., C seg and C cr ) in an adversarial scheme to reduce domain shifts at local scales to achieve local contextual-relation consistency.

Contextual-relation consistent domain adaptation
This subsection introduces our contextual-relation consistent domain adaptation at local scales, denoted as CrCDA * , via adaptive entropy max-minimizing, as shown in Fig. 2.
Contextual-relation pseudo label establishment. In order to implement local-scale task, we sample regions on the feature space and implement domain alignment at local scales to achieve local contextual-relation consistent domain adaptation, as shown in Fig. 3. Different from [22] that implements mode-agnostic patches alignment or [45] that aligns patch-indexed representation of the whole image only at global scales by a discriminator (i.e., the probability distributions of patch index prediction of the whole images.), we aim to aligns inter-class relations within each single patch, i.e., the probability distributions of pixel class prediction within each patch, w.r.t its mode via a classifier. Thus, the preliminary is to establish the region-scale label, where we first crop the pixel-scale ground-truth to many larger regions and then use Dbscan [12] to cluster them to assign each region a certain index label (i.e., contextual-relation pseudo label). Specifically, we assign the index label to regions according to the clustering results based on the histogram of gradient. For region-scale label (i.e., contextual-relation pseudo label), we cluster regions into different groups based on the histogram of gradient and assign the index label. These region-scale/contextual-relation pseudo labels can assist our network to implement alignment at local scales. Detailed information about the regionscale/contextual-relation pseudo labels is in the supplementary materials.  Adaptive entropy max-minimizing adversarial learning scheme. In local-scale adaptation, C cr aims to approximate the prototypical feature representations for each contextual-relation (e.g., road-sidewalk, sky-building, polesidewalk, etc.) by implementing entropy maximization in target domain according to the source prototypical feature representations found via supervised learning in source domain. E focuses on extracting discriminative feature representations (near the approximated prototypical feature representations) by implementing entropy minimization. Specifically, the prototypical feature representations of source domain found with supervision are first utilized to estimate the prototypical feature representations for target data by entropy maximizing w.r.t C cr . E then adapts the extracted feature representations to the corresponding prototypical feature representations by minimizing the entropy. The overall unsupervised domain adaptation at local scales is achieved by the adversarial training between C cr and E as illustrated in Fig. 4. Different from that applied in semi-supervised learning [37], our unsupervised domain adaptation training method, referred as adaptive entropy max-minimizing (AEMM) implements entropy max-min with a regularizer R(P ) = ave{P log P } × λ R (λ R decreases with training iteration, details are shown in appendix) for better estimating the prototypes in the target domain where no labels are available.
Source Flow. In our local-scale adaptation setting, the source data contributes to L seg and L cr . Given a source image x s ⊂ X s , its corresponding segmentation label y s ⊂ Y s and contextual-relation pseudo-label y s cr ⊂ Y s cr , P (h,w,c) s = C seg (E(x s )) is the predicted probability map w.r.t each pixel over C classes; P (i,j,n) s cr = C cr (E(x s )) is the predicted probability map w.r.t each region over N pre-defined contextual-relations classes. Therefore, it is a simple supervised learning objective to minimize L seg and L cr , which are expressed as: Target Flow. As the target label is not accessible, we introduce the adversarial training scheme between feature extractor E and classifier C cr to extract discriminative features for target data via adaptively max-minimizing entropy in target domain. Given a target image x t ⊂ X t , P (i,j,n) t cr = C cr (E(x t )) is the predicted probability map w.r.t each region over N pre-defined contextual-relations classes. The entropy loss L entcr is expressed as: For local-scale adaptation, we use the same back-propagation optimizing scheme with the gradient reverse layer mentioned in [57]. The training objective can be express as: where λ ent is a weight factor to control the balance of unsupervised adaptation on target data and supervised learning on source data.

CrCDA with pixel-/global-scale
This subsection introduces our CrCDA with pixel-/global-scale, denoted as Cr-CDA, via adaptive entropy max-minimizing, as shown in Fig. 2. Our discriminatorfree AEMM adversarial training scheme can also be extended into pixel-scale and global/image-scale to form multi-scale domain adaptation.
In multi-scale adaptation, for L seg , L cr and L ent cr , the objectives are the same as that in local-scale adaptation. We extend the AEMM adversarial training scheme mentioned before into pixel-scale and global-scale adaptation. For pixel-scale adaptation, we implement pixel-scale entropy loss L ent on target data to E and C seg . For global-scale adaptation, we implement global-scale entropy loss L ent D on target data to E and C D , where C D is a domain classifier. C D takes the layout probability map concatenated by the two probability maps generated from C seg and C cr as input, and predicts domain label for it (e.g., 0 for source domain, 1 for target domain). The global-alignment is achieved by the adversarial training between C D and (E, C seg , C cr ). Finally, our multi-scale consistent domain adaptation network is able to align domain shift at global scales, local-scale and pixel-scale simultaneously.
Similar to local-scale adaptation, we formulate the pixel-scale entropy loss as: t pix ), 0} (5) For multi-scale adaptation, we also use the same back-propagation optimizing scheme with the gradient reverse layer mentioned in [13,14]. The training objective can be express as: where L D is provided in supplementary materials; λ cr , λ ent and λ D are the weight factor to balance the unsupervised adaptation on target data and the task-specific objectives on source data.

Datasets
We evaluate our unsupervised domain adaptation networks for semantic segmentation on two challenging synthesized-to-real tasks: GTA5 [35] → Cityscapes [9] and SYNTHIA [36] → Cityscapes. GTA5 contains 24, 966 synthesized images with high-resolution and we use the 19 common categories between GTA5 and Cityscapes in the same setting as in [44]. SYNTHIA contains 9, 400 synthetic images with 16 common categories in Cityscapes. We use either GTA5 or SYN-THIA as source domain. We use the unlabelled training set of Cityscapes as target domain, which includes 2975 real-world images.

Implementation Details
For a fair comparison, similar to [49] [31] [44], we utilize Deeplab-V2 architecture [3] with ResNet-101 pretrained on ImageNet [10] as our single-scale semantic Table 1. Results of domain adaptation task GTA5 → Cityscapes. "V" means the VGG16-based model and "R" means the ResNet101-based model.  segmentation network (E + C seg ). To extend our model to multi-scale network, we simply copy and modify C seg to create C cr and C D with different output channels (e.g.,N and 1) and different output sizes due to various scales (i.e., region-size and global-size). We also apply our methods on VGG-16 [43] in the same way as employing ResNet-101. Following [13] [47], a gradient reverse layer is employed to reverse the entropy loss between E and (C seg , C cr ) during pixel-/region-scale adaptation to achieve adversarial training. The domain classifier C D works similar to a discriminator for global-scale alignment. During training, we utilize SGD [1] to optimize our networks with a momentum of 0.9 and a weight decay of 1e − 4. The initial learning rate is set as 2.5e − 4 and decayed by a polynomial policy with a power of 0.9, as illustrated in [3]. For all experiments, the hyper-parameters λ ent , λ D , λ cr and N are set as 2.5e − 5, 2.5e − 5, 5e − 3 and 100, respectively.

Comparison with state-of-art
We compare the experimental results of our method and state-of-the-art algorithms in two "Synthetic-to-real" UDA tasks with two different architectures: VGG-16 and ResNet-101. For "GTA5 → Cityscapes", we present the results in Table 1 with comparisons to the state-of-the-art domain adaptation methods [49,31,44,19,55,18,56,57]. Our contextual-relation consistent domain adaptation, expressed as CrCDA, achieves comparable performance to other stateof-the-art approaches on both architectures. Compared to Adapt-SegMap (output space global alignment) [44], category-level adversarial network (output space class-wise alignment) [31] and patch-represented global alignment [45] (patch-indexed latent space alignment), CrCDA consistently brings over +2.1% mIoU improvements on ResNet-101. We reckon this gain is from our end-toend/concurrent multi-scale alignment, which indicates that local consistency Table 2. Results of domain adaptation task SYNTHIA → Cityscapes. "V" means the VGG16-based model and "R" means the ResNet101-based model. "mIoU" and "mIoU*" are calculated over 16 and 13 classes, respectively.  Table 2, we present the adaptation result for the task "SYNTHIA → Cityscapes" and consistent improvements are observed w.r.t state-of-the-arts. Detailed analysis is included in next subsection.

Ablation Studies and Analysis
We analyze our proposed CrCDA with several state-of-the-art baselines. In general, both single-scale form (CrCDA * ) and multi-scale form (CrCDA) achieve comparable results to all the baselines in all the settings.
As shown on the first three rows in Table 3, our pixel-scale AEMM adversarial network brings +1.4% improvements in terms of mIoU over MinEnt [49]. The reason lies in that direct entropy minimization does not take the domain gap into account while our AEMM training scheme pushes the source distribution closer to target distribution during maximizing entropy on target data.
For our CrCDA with single-scale form (CrCDA * ) via AEMM, it outperforms MinEnt-based contextual-relations alignment by +1.6% on ResNet-101, as shown on the second block (row4-5) in Table 3. We reckon that these improvements are contributed by our adaptive entropy max-min training scheme which considers the domain mismatch/gap while MinEnt neglects.
Our CrCDA with multi-scale form integrating three scales' adaptation (pixel-, local-and global-scale), termed as CrCDA shown on the bottom block in Table 3, achieves state-of-the-art performances 48.6% mIoU on ResNet-101. Besides, CrCDA also outperforms all current methods by over +1.5%. Compared to "Pixel+Global", CrCDA brings +2.6% improvement in mIoU, which demonstrates that local-scale alignment is essential as well as other scales (e.g., pixelscale and global-scale). In fact, the local contextual-relation consistent adaptation loss (i.e., L ent cr ) penalize groups of pixels predictions to achieve local-scale Table 3. Ablation study of the proposed contextual-relation consistent domain adaptation on GTA5-to-Cityscapes using the ResNet-101 network. All settings/methods are with "Lseg" (bold texts represent our methods). CrCDA * represents the contextualrelation consistent domain adaptation with only single-scale (local scale).   [45] 43.2 CrCDA * + AdaptSeg [44] 44.8 CrCDA * + PatAlign [45] 44.7 CrCDA * + AdvEnt [49] 45.2 CrCDA * +GlobalAlign (ours) 46.1 alignment, where global-scale adaptation loss operates more on image-scale (e.g., scene layout) while that of pixel-scale works on the feature representation alignment of each independent pixels. The consistent results with different settings further confirm that complementary information has been learned in different scales' adaptation. The qualitative results and visualization of feature distributions are provided in Fig. 5 and 6, which further demonstrate our conjectures mentioned above. We also provide the complementary studies to demonstrate that our local contextual-relations alignment method is complementary to most existing global-scale alignment approaches, as shown in Table 4.

Conclusions
In this paper, we present the local contextual-relation consistent domain adaptation (CrCDA) to address the task of unsupervised domain adaptation for seman-  Fig. 6. Visualization of feature distributions via t-SNE [33]. "Blue": Source. "Red": Target. As shown in the first column, the feature distribution of source data is naturally more discriminative (discrete) than that of target data (uniformly distributed) due to only source supervision is available. Traditional global alignment (TGA) aligns them in global scale, where global consistency is achieved while local consistency is ignored. Thus the adapted target feature distribution is not discriminative. CrCDA aligns them with local-scale consistency (i.e., local contextual-relation consistency), where both local and global consistencies are achieved. Thus the adapted target feature distribution is more discriminative and consistent with that of the source.
tic segmentation. By taking a closer look at the local inconsistency (i.e., local contextual-relations inconsistency) while implementing global adaptation, Cr-CDA is able to align the domain shift in local and global scales at the same time, where local semantic consistency is normally ignored by current approaches. The experimental results on the two challenging segmentation UDA tasks validate the state-of-the-art of CrCDA.

Acknowledgement
This research was conducted in collaboration with Singapore Telecommunications Limited and partially supported by the Singapore Government through the Industry Alignment Fund -Industry Collaboration Projects Grant.

A. Appendix
A.1. Method -Details of region-scale/contextual-relation pseudo labels and regularizer weight We would share more details about the region-scale/contextual-relation pseudo labels and the weight of regularizer used in this paper. For the source domain, the sizes of the input image for datasets GTA5 and SYNTHIA are 720 × 1280 and 760 × 1280 , respectively. In this paper, we use two types of regions with two different sizes. The first sizes of regions for datasets GTA5 and SYNTHIA are 18 × 32 and 19 × 32, respectively. The second sizes of regions for datasets GTA5 and SYNTHIA are 36 × 64 and 38 × 64, respectively. For the target domain (dataset Cityscapes), the size of input image is 512 × 1024. The sizes of regions are 16×32 and 32×64, respectively. We use two independent contextual-relations (CR) classifiers to deal with these two types of regions with two different sizes.
The weight of the regularizer in adaptive entropy max-minimizing adversarial learning scheme decreases with training iteration, which is expressed as: λ R = (1 − iter max iter ) power with power = 0.9.

A.2. Method -Traditional Losses
For the source domain, traditional approaches learn a supervised segmentation model G that aims to minimize a segmentation loss. For the target domain, UDA networks using adversarial learning train G to extract domain-invariant features though the minimaxing game between G and a domain discriminator D. The overall loss in the UDA networks can therefore be formulated by: L(X s , X t ) = L seg (G) + L adv (G, D) A.3. Method -Loss in Multi-Scale Adaptation Source Flow: In our contextual-relation consistent domain adaptation (Cr-CDA) with multi-scale form, the source-domain data contribute to L seg , L cr and L D . Given a source-domain image x s ⊂ X s and the corresponding pixelscale label y s ⊂ Y s , region-scale (contextual-relations) pseudo label y s cr ⊂ Y s cr , P (h,w,c) s = C seg (E(x s )) is the predicted probability map w.r.t each pixel over C classes; P (i,j,n) s cr = C cr (E(x s )) is the predicted probability map w.r.t each region over N classes. The layout probability map P (h,w,c+n) s layout is generated by concatenating P Target Flow: As the target label is not accessible, we design an adversarial training scheme between feature extractor E and classifiers (C seg , C cr and C D ) that extracts discriminative features via max-minimizing entropy in the target domain. Given a target image x t ⊂ X t , P (h,w,c) t = C seg (E(x t )) is the predicted probability map w.r.t each target pixel over C classes; P (i,j,n) t cr = C cr (E(x t )) is the predicted probability map w.r.t each target region over N classes. The layout probability map P Therefore, the overall global alignment loss is expressed as: where domain classifier entropy is Ent s d = −C D (P

A.4. Experiment -More Qualitative Results
We share more qualitative experimental results for GTA5 → Cityscapes as shown in Fig. 7. As Fig. 7 shows, our CrCDA aligns both low-level features (e.g., boundaries of sidewalk, car and person etc.) and high-level features by multi-scale adversarial learning. As a comparison, AdvEnt neglects low-level information which focuses more on high-level features. As a result, CrCDA achieves both local and global consistencies in segmentation while AdvEnt achieves global consistency only. Qualitative results for GTA5 → Cityscapes. Our approach (CrCDA) aligns low-level features (e.g., boundaries of sidewalk, car and person etc.) as well as highlevel features by multi-scale adversarial learning. In contrast, AdvEnt ignores low-level information because global alignment focuses more on high-level information. Thus, as shown above, CrCDA achieves both local and global consistencies while AdvEnt only achieves global consistency.