Nikola Banovic
email: nbanovic@cs.cmu.edu

Tofi Buzali

Fanny Chevalier
email: fanny.chevalier@inria.fr

Jennifer Mankoff
email: jmankoff@cs.cmu.edu

Anind K Dey

Modeling and Understanding Human Routine Behavior

Keywords: Inverse Reinforcement Learning; Markov Decision Process, ACM Classification H, 5, m, Information interfaces and presentation (e, g, , HCI): Miscellaneous

Human routines are blueprints of behavior, which allow people to accomplish purposeful repetitive tasks at many levels, ranging from the structure of their day to how they drive through an intersection. People express their routines through actions that they perform in the particular situations that triggered those actions. An ability to model routines and understand the situations in which they are likely to occur could allow technology to help people improve their bad habits, inexpert behavior, and other suboptimal routines. However, existing routine models do not capture the causal relationships between situations and actions that describe routines. Our main contribution is the insight that byproducts of an existing activity prediction algorithm can be used to model those causal relationships in routines. We apply this algorithm on two example datasets, and show that the modeled routines are meaningful-that they are predictive of people's actions and that the modeled causal relationships provide insights about the routines that match findings from previous research. Our approach offers a generalizable solution to model and reason about routines.

INTRODUCTION

Routine behavior defines the structure of and influences almost every aspect of people's lives. Routines are defined by frequent actions people perform in different situations that are the cause of those actions [START_REF] Hodgson | The ubiquity of habits and rules[END_REF]. Routines are a type of purposeful behavior made up of goal-directed actions [START_REF] Taylor | Purposeful and non-purposeful behavior: A rejoinder[END_REF], which people acquire, learn and develop through repeated practice [START_REF] Martin | Habit: the 95% of behavior marketers ignore[END_REF][START_REF] Ronis | Attitudes, decisions, and habits as determinants of repeated behavior[END_REF]. As such, good routines enable predictable and efficient completion of frequent and repetitive tasks and activities [START_REF] Orlikowski | It's About Time: Temporal Structuring in Organizations[END_REF]. Routines can describe people's daily commute, their sleeping and exercising patterns, or even low-level tasks, such as how they operate their vehicle through an intersection. Routines, like most other kinds of human behaviors, are not fixed, but instead may vary and adapt based on feedback and preference [START_REF] Feldman | Reconceptualizing organizational routines as a source of flexibility and change[END_REF].

A key aspect of being able to understand and reason about routines is being able to model the causal relationships between people's situations and the actions that describe the routines. We refer to frequent departures from established routines as routine variations, which are different from deviations and other uncharacteristic behavior that do not contribute to the routines. An ability to model routines and their variations, within and across people, could help researchers better understand routine behavior, and inform technology that influences routine behavior and helps people improve the quality of their lives.

Studies of routines often characterize routines in terms of a series of actions, typically derived from large activity data sets. Data mining algorithms may automatically extract patterns from such data (e.g., [START_REF] Brdiczka | Temporal task footprinting: identifying routine tasks by their temporal patterns[END_REF][START_REF] Eagle | Eigenbehaviors: identifying structure in routine[END_REF][START_REF] Farrahi | Extracting Mobile Behavioral Patterns with the Distant N-Gram Topic Model[END_REF][START_REF] Magnusson | Discovering hidden time patterns in behavior: T-patterns and their detection[END_REF][START_REF] Sadilek | Far Out: Predicting Long-Term Human Mobility[END_REF]), while visualization can help researchers to interrogate the data (e.g., [START_REF] Aigner | PlanningLines: Novel Glyphs for Representing Temporal Uncertainties and Their Evaluation[END_REF][START_REF] Plaisant | LifeLines: visualizing personal histories[END_REF][START_REF] Wongsuphasawat | LifeFlow: visualizing an overview of event sequences[END_REF]). However, these existing approaches do not model the causal relationship between situations and actions. This makes it difficult to study and explain which situations or contexts, defined as the environmental information relevant to an individual's current activity [START_REF] Anind | Understanding and Using Context[END_REF], trigger which routine actions. A different approach is to detect [START_REF] Bulling | A tutorial on human activity recognition using bodyworn inertial sensors[END_REF], classify [START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF] and predict [START_REF] Baratchi | A hierarchical hidden semi-Markov model for modeling mobility data[END_REF][START_REF] Krumm | Predestination: inferring destinations from partial trajectories[END_REF] activities from large data sets. Although such approaches imply the causality between the contexts and actions, it is difficult to extract meaning about routines from models learned using those algorithms and understand the reasons why they make certain classifications and predictions [START_REF] Lim | Investigating intelligibility for uncertain context-aware applications[END_REF].

To address those issues, we present a novel approach to automatically extract and model routines and routine variations from human behavior logs. Our approach supports both individual and population models of routines, providing the ability to identify the differences in routine behavior across different people and populations. Our main contribution to modeling routines is our insight that the byproducts of MaxCausalEnt [START_REF] Ziebart | Modeling purposeful adaptive behavior with the principle of maximum causal entropy[END_REF], a decision-theoretic algorithm typically used to predict people's activity, actually encode the causal relationship between routine actions and context in which people perform those actions. These causal relationships allow for reasoning about and understanding of the extracted routines.

Using two different existing human activity data sets, we evaluate the ability of our approach to extract different types of routines from diverse types of behavior: people's daily schedules and commutes [START_REF] Davidoff | How routine learners can support family coordination[END_REF] and activities that describe how people operate a vehicle [START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF]. We show that the extracted routine patterns are at least as predictive of behaviors in the two behavior logs as the baseline we establish with existing algorithms. Next, we recruited researchers that work with human activity and routine data to verify that patterns extracted using our approach are meaningful and match the ground truth reported in previous work [START_REF] Davidoff | How routine learners can support family coordination[END_REF][START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF]. For the purposes of this task, we developed a tool that enabled the researchers to visually explore and compare the routines extracted using our approach.

Our first set of results show that our approach enables extraction of a reasonable set of human readable patterns of routine behavior from behavior logs. Our second set of results demonstrate that the models of causal relationships provided by our approach can help researchers explore, understand, and form new insights about human routines from behavior logs without having to manually search for those patterns in raw data. Another benefit of our approach compared to existing systems is that it promises efficient automated prediction and reasoning about routines even under uncertainty that is inherent in human behavior. We discuss how such models and knowledge about routines can inform the design of novel systems that help people improve their behavior.

UNDERSTANDING HUMAN ROUTINE BEHAVIOR

Many different stakeholders may care about understanding routine behavior. For example, a clear picture of the many aspects of routine behavior can help researchers to generate theories and models of human routine behavior. Designers may build on such theories to design technologies that help people to improve their routines [START_REF] Davidoff | How routine learners can support family coordination[END_REF][START_REF] Davidoff | Learning patterns of pick-ups and drop-offs to support busy family coordination[END_REF]. Models of human routines may be used for prediction and automation [START_REF] Ziebart | Navigate like a cabbie: probabilistic reasoning from observed context-aware behavior[END_REF]. Individuals may wish to reflect about their own routines for better understanding and supporting behavior change [START_REF] Li | A stagebased model of personal informatics systems[END_REF].

Although routine behaviors result from low-level cognitive plans, which can be modeled using existing cognitive architectures (e.g., [START_REF] Anderson | An integrated theory of the mind[END_REF]), such models are difficult to apply in understanding high-level routine behaviors. Thus our focus in this work is on how routines are expressed in actions people perform in different contexts.

Through an analysis of stakeholder use cases, and a review of the literature, we derive three research questions that these stakeholders, and in particular researchers, who analyze routines may wish to answer in order to meet their goals. We review existing approaches to understanding routines based on how well they answer those questions.

Current Approaches to Understanding Routine Behavior

Visualizing data from behavior logs is a common way for researchers to identify routines. Logged behavior data is often visualized on a timeline as a sequence of events. The simplicity of this approach makes it applicable to a variety of domains, such as schedule planning to show uncertainty of duration of different events [START_REF] Aigner | PlanningLines: Novel Glyphs for Representing Temporal Uncertainties and Their Evaluation[END_REF], visualizing family schedules [START_REF] Davidoff | How routine learners can support family coordination[END_REF], and representing events related to patient treatments [START_REF] Plaisant | LifeLines: visualizing personal histories[END_REF]. More advanced timelines enable the user to specify properties of the timeline for easier viewing. For example, Spiral Graph [START_REF] Weber | Visualizing time-series on spirals[END_REF] aligns sequential events on a spiral timeline using a user-defined period. However, due to the complexity and size of behavior logs, simply visualizing raw event data does not guarantee that the user will be able to find patterns of behaviors that form routines.

Other visualization approaches enable users to find patterns in time-series data by manually querying and highlighting different parts of behavior event sequences [START_REF] Buono | Interactive pattern search in time series[END_REF][START_REF] Zhao | KronoMiner: using multi-foci navigation for the visual exploration of time-series data[END_REF][START_REF] Zhao | Exploratory analysis of time-series with chronolenses[END_REF] or manually aggregating common sequences of events based on features in the data [START_REF] Jing | Interactive querying of temporal data using a comic strip metaphor[END_REF][START_REF] Monroe | Temporal Event Sequence Simplification[END_REF] until meaningful patterns emerge. The user is then able to judge the quality and saliency of the patterns by visual inspection. However, during the early exploratory stages, users might not always know what features are important and contribute to routine behavior patterns, making manual exploration challenging.

One major limitation of existing visualizations lies in their lack of support for both context and actions. Rather, they focus on isolated events, or temporal evolution of a particular state (e.g., sleep vs. awake) or variable (e.g., the amount of steps walked per day). Visualizing both context and actions is challenging partially because even advanced interactive visualizations have difficulty in visualizing routine patterns that depend on multiple heterogeneous variables, especially as the number of variables grows. Automated routine extraction and summarization is another option for exploring routines in behavior logs. However, patterns extracted using the existing methods often do not include important aspects of routine behavior. For example, T-patterns [START_REF] Brdiczka | Temporal task footprinting: identifying routine tasks by their temporal patterns[END_REF][START_REF] Magnusson | Discovering hidden time patterns in behavior: T-patterns and their detection[END_REF] can automatically find recurrences of events in behavior logs, but do not explicitly capture the contexts and actions. Other methods, such as Topic Models [START_REF] Farrahi | Extracting Mobile Behavioral Patterns with the Distant N-Gram Topic Model[END_REF], include features that describe both context and actions, but without modeling the structure of possible variations from those routines. Methods based on Hierarchical Task Networks [START_REF] Li | Learning Probabilistic Hierarchical Task Networks as Probabilistic Context-Free Grammars to Capture User Preferences[END_REF] and Eigen decomposition [START_REF] Eagle | Eigenbehaviors: identifying structure in routine[END_REF][START_REF] Sadilek | Far Out: Predicting Long-Term Human Mobility[END_REF] capture the structural components of the contexts and actions, but do not explicitly model the causal relationship between the two that defines the routine [START_REF] Hodgson | The ubiquity of habits and rules[END_REF]. While these algorithms are helpful for extracting routine, they are not sufficient for helping to understand routine behaviors. This leads to our first research question that researchers want to answer: RQPat: How can we expose the full complexity of routine behavior patterns? To extract and understand routines, it is critical to discern which features of the context influence routine actions and which features of the actions people have a demonstrated preference for.

Where Routines End and Variations Begin

How people respond to routine variations is also an important part of routine behavior. Human routine behavior is not static, which makes variations from the routines inevitable [START_REF] Feldman | Reconceptualizing organizational routines as a source of flexibility and change[END_REF]. Variations often result from new contexts or unforeseen circumstances for which people have not yet found a routine [START_REF] Davidoff | How routine learners can support family coordination[END_REF], or occasions when people want to explore alternative ways to accomplish their tasks [START_REF] Feldman | Reconceptualizing organizational routines as a source of flexibility and change[END_REF]. For example, parents may vary their routine when their child has a new scheduled activity that is not part of their current routines or when they unexpectedly have to pickup their child from an existing activity [START_REF] Davidoff | How routine learners can support family coordination[END_REF]. Understanding variations is also important to understand the tradeoffs between different behaviors that people have. For example, in designing systems that help people comply with their gym routine, it is important to understand what other activities cause people to depart from their gym routine.

However, existing machine learning algorithms (e.g., [START_REF] Baratchi | A hierarchical hidden semi-Markov model for modeling mobility data[END_REF][START_REF] Bulling | A tutorial on human activity recognition using bodyworn inertial sensors[END_REF][START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF][START_REF] Krumm | Predestination: inferring destinations from partial trajectories[END_REF]) purposefully disregard variations in human behavior to focus on classifying and predicting only the most frequent human activity. Also, some variations may happen infrequently in data and are difficult to detect using those existing algorithms. Some specific infrequent variations may be detectable (such as detecting when parents are going to be late to pickup their children [START_REF] Davidoff | Learning patterns of pick-ups and drop-offs to support busy family coordination[END_REF]). However, this requires a case-by-case approach to address each kind of variation, which can be difficult to apply if all possible variations are not known a priori. This leads to our second research question researchers want to answer: RQVar: How can we expose variations of routines and enable the discovery of the cause of those behaviors? Similar to our first question, routine-modeling approaches must identify aspects of context that influence how people behave in the face of variations.

Routines Across Individuals and Populations

The third issue to consider is how to explore the differences in routines within and across individuals and populations. Different people often develop their own routines to deal with their individual contexts. For example, differences in family daily routines can explain the responsibilities of different family members [START_REF] Davidoff | How routine learners can support family coordination[END_REF]. Also, to help people improve their routines it is important to understand the difference between peoples' desired routines and their actual routines. For example, to design systems that help people better organize their schedules to avoid being late to pickup their children requires understanding differences between routines for days when they are late and days when they are on time.

Such comparisons can be performed between routines of an individual, but also between routines of different people or populations. For example, to create interventions that help aggressive drivers improve their driving style requires an understanding of the differences between aggressive and non-aggressive driving routines [START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF]. Based on how close drivers are to aggressive or non-aggressive routines, they can be classified into aggressive and non-aggressive groups. Researchers can then identify routine behaviors that need to change to make aggressive drivers less aggressive overall.

Existing algorithms can, of course, be applied to individual or population models of routines. For example, a standard strategy is to build a model for each population and then compare them to establish the differences and similarities between the populations. Population models can then be visualized using existing approaches: Hierarchical Task Networks [START_REF] Li | Learning Probabilistic Hierarchical Task Networks as Probabilistic Context-Free Grammars to Capture User Preferences[END_REF] can be visualized using existing Probabilistic Context-Free Grammar-based tools (e.g., [START_REF] Collins | Visualization of uncertainty in lattices to support decision-making[END_REF]), and Arc Diagrams [START_REF] Wattenberg | Arc diagrams: Visualizing structure in strings[END_REF] can visualize temporal patterns extracted using T-patterns [START_REF] Magnusson | Discovering hidden time patterns in behavior: T-patterns and their detection[END_REF]. Those population models can then be coordinated and displayed across multiple views [START_REF] Roberts | State of the Art: Coordinated & Multiple Views in Exploratory Visualization[END_REF]. However, such visualization techniques are not yet widely adopted. Finding differences in models is still often based on intuition and expert domain knowledge, and limited to visually comparing the distributions of feature values. This leads to our third question researchers want to answer: RQComp: How can we support comparison of routines across individuals and populations? Therefore, we need to ensure that different aspects of routines and routine variations that support comparison between different people and populations are included in the models.

MODELING HUMAN ROUTINE BEHAVIOR

In this section we present our approach to modeling human routine behavior that applies an existing decision-theoretic algorithm to the domain of routine modeling. To identify routine patterns (RQPat), we explicitly model the causal relationship between contexts in which different routines occur and the actions that people perform in those contexts. Unlike models that extract only the most frequent routines, our approach also models possible variations from those routines (RQVar), even in infrequent contexts. Our approach does this by modeling probability distributions over different possible behaviors, which allows the researcher to make sense about which of those behaviors form routines and which form variations. Our approach can model both individual and population routine behavior, and thus allows comparisons between those models (RQComp).

Data Modeling

Human behavior data is often collected using different sensors and stored into behavior logs. In our approach we first convert the behavior logs into sequences of events representing people's current context and the actions they perform in that context. We then use those sequences of events to model human routine behavior.

Model of Human Routine Behavior

We model demonstrated routine behavior using a Markov Decision Processes (MDP) framework [START_REF] Puterman | Markov decision processes: discrete stochastic dynamic programming[END_REF]. MDP is particularly well suited for modeling human routine behavior because it explicitly models the user's context, the actions that can be performed in that context, and the preferences people have for different actions in different contexts. A Markov decision process is a tuple: It consists of a set of states representing context, and actions that a person can take. In addition, the model includes an action-dependent probability distribution for each state transition , which specifies the probability of the next state when the person performs action in state . This state transition probability distribution models how the environment responds to the actions that people perform in different states. When modeling human behavior, the transitions are often stochastic (each pair can transition to many transition states with different probabilities). However, if the person has full control over the environment, they can also be deterministic (i.e., for each pair there is exactly one transition state with probability 1). Finally, there is a reward function that the person incurs when performing action in state , which represents the utility that people get from performing different actions in different contexts.

People's behavior is then defined by sequences of actions they perform as they go from state to state until reaching some goal state. In an MDP framework, such behavior is defined by a deterministic policy ( ), which specifies actions people take in different states. Traditionally, the MDP is "solved" using algorithms, such as value iteration [START_REF] Bellman | A Markovian decision process[END_REF], to find an optimal policy (with the highest expected cumulative reward). However, our goal is to find the expected frequencies of different states and the probability distribution of actions given states instead-information necessary to identify people's routines and variations.

Learning Routine Patterns from Demonstrated Behavior

In this section we explore how the MaxCausalEnt algorithm [START_REF] Ziebart | Modeling purposeful adaptive behavior with the principle of maximum causal entropy[END_REF], an algorithm typically used to predict human behavior [START_REF] Ziebart | Navigate like a cabbie: probabilistic reasoning from observed context-aware behavior[END_REF][START_REF] Ziebart | Planning-based prediction for pedestrians[END_REF], can be applied in a novel way to extract routine behavior patterns from observed data. MaxCausalEnt algorithm makes its predictions by computing a policy ( ) that best predicts the action people take in different states. Our main contribution is our insight that in the process of computing this policy, MaxCausalEnt algorithm computes two other functions that express how likely it is that a state and action are part of a routine: 1) the expected frequency of states ( ), and 2) probability distribution of actions given states . We now describe how we compute these two functions and how they relate to routines.

Inverse Reinforcement Learning (IRL) [START_REF] Ng | Algorithms for Inverse Reinforcement Learning[END_REF] approaches, which MaxCausalEnt is based on, assume that people assign a utility function (modeled as the reward functions

), which they use to decide which action to perform in different demonstrated states. Each state and action combination in our MDP model is expressed by a feature vector . For example, in an MDP that models daily commute routines, states can have features that describe all possible locations that a person can be at, and actions can have features that describe if the person is staying at or leaving the current location. As is common for IRL algorithms [START_REF] Ng | Algorithms for Inverse Reinforcement Learning[END_REF][START_REF] Ziebart | Modeling interaction via the principle of maximum causal entropy[END_REF], we assume a parametric reward function that is linear in , given unknown weight parameters :

We begin the process of recovering the expected state frequencies ( ) and probability distribution of actions given states (

) by trying to learn the person's reward functions from demonstrated behavior. This problem reduces to matching the model feature function expectations (

) with demonstrated feature expectations ( [START_REF] Abbeel | Apprenticeship learning via inverse reinforcement learning[END_REF]. To match the expected counts of different features, we use MaxCausalEnt IRL [START_REF] Ziebart | Modeling interaction via the principle of maximum causal entropy[END_REF], which learns the parameters of the MDP model to match the actual behavior of the person. Unlike other approaches described earlier, MaxCausalEnt explicitly models the causal relationships between context and actions, and keeps track of the probability distribution of different actions that people can perform in those contexts.

To compute the unknown parameters , MaxCausalEnt considers the causal relationships between all the different features of the states and the actions. The Markovian property of MDP, which assumes that the actions a person performs only depend on the information encoded by the previous state, makes computing the causal relationships between the states and actions computationally feasible. MaxCausalEnt extends the Principle of Maximum Entropy [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] to cases where information about probability distribution is sequentially revealed, as is the case with behavior logs. This principle ensures that the estimated probability distribution of actions given states is the one that best fits the state and action combinations from the sequences in the behavior logs.

MaxCausalEnt IRL maximizes the causal entropy ( ) of the probability distribution of actions given states :

such that:

The first constraint in the above equation ensures that the feature counts calculated using the estimated probability distribution of actions given states ( ) matches the observed counts of features in the data, and the other two ensure that is an actual probability distribution.

Using the action-based cost-to-go ( ), which represents the expected value of performing action in state , and state-based value ( ) notation, which represents the expected value of being in state , the procedure for MDP MaxCausalEnt IRL reduces to [START_REF] Ziebart | Modeling interaction via the principle of maximum causal entropy[END_REF]: Note that this is similar, but not the same as stochastic value iteration [START_REF] Bellman | A Markovian decision process[END_REF], which would model optimal and not observed behavior. The probability distribution of actions given the states is then given by:

The probability distribution of actions given states and the state transition probability distribution are used in a forward pass to calculate the expected state frequencies ( ). This optimization problem can then be solved using a gradient ascent algorithm. Ziebart [START_REF] Ziebart | Modeling purposeful adaptive behavior with the principle of maximum causal entropy[END_REF] provides proofs of these claims and detailed pseudocode for the algorithm above.

Extracting Models of Routine Behavior

We illustrate our routine modeling approach on two previously collected data sets from the literature that contain logs of demonstrated human behavior. The first data set contains daily commute routines of all family members from three two-parent families with children from a mid-sized city in North America [START_REF] Davidoff | How routine learners can support family coordination[END_REF]. The data set was used to predict the times the parents are likely to forget to pickup their children [START_REF] Davidoff | Learning patterns of pick-ups and drop-offs to support busy family coordination[END_REF]. The other data set contains driving routine behavior of aggressive and non-aggressive drivers as they drive on their daily routes [START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF]. The data set was used to classify aggressive and non-aggressive drivers.

We picked these two data sets to show the generalizability of our approach to different types of routines. The two data sets contain routine tasks people perform on a daily basis, but that are very different in nature. The family daily routine data set incorporates the traditional spatio-temporal aspect of routines most of the existing work focuses on. The driving data set contains situational routines that are driven by other types of context (e.g., the surrounding traffic, the current position of the car in the intersection).

The two data sets also differ in granularity of the tasks. The commute routines happen over a longer period of time and the granularity of the task is very coarse with few actions that people can perform in different contexts (e.g., stay at the current place or leave and go to another place). The daily routines are therefore defined by the states the people are in. The aggressive driving data set contains fine-grained actions, which often occur in parallel, that people perform to control the vehicle (e.g., control the gas and brake pedals and the steering wheel). Driving routines are therefore primarily defined by the drivers' actions in different driving situations. The driving data set also showcases the ability of our approach to capture population models (e.g., aggressive drivers vs. non-aggressive drivers) and enable comparison of routines across different populations.

The data in our two example data sets consists of sequences of sensor readings, which we convert into sequences of events represented by state action pairs. Parsing the raw data we are able to extract: 1) a set of states defined by a set of features which represent context, 2) a set of actions defined by a list of binary features which represent activities that the people can perform, and 3) empirically estimated state-action transition dynamics (

). At any discrete event step, the state features contain values of all the contextual sensor readings at that event, and actions contain feature values describing the activity the people performed at that event. We then estimate state-action transition dynamics based on the frequencies of state transitions in the state-action event sequences and estimate the expected state frequency counts ( ) and the state-action probability distributions ( ) as described in the previous section.

Family Daily Routines Data Set

Situations when one of the parents is unable to pickup or drop-off a child create stress for both parents and children [START_REF] Davidoff | How routine learners can support family coordination[END_REF]. To better understand the circumstances under which these situations arise, it is important to identify when the parents are responsible for picking up and dropping off their children (RQPat), when variations from normal routines occur and how parents handle those situations (RQVar). This requires finding and understanding how the parents organize their daily routines around those pickups and drop-offs (RQComp). This data contains location sampling (latitude and longitude) at one-minute intervals for every family member (including children) in three families from a mid-sized city in North America [START_REF] Davidoff | How routine learners can support family coordination[END_REF]. Location information was manually labeled based on information from bi-weekly interviews with participants. Participants also provided information about their actual daily routines during those interviews.

We converted the location logs into sequences of states and actions representing each individual's daily commute for each day in the data set. State features included the day of the week, hour of the day, participant's current place, and whether the participant stayed at the location from the previous hour, arrived at the location during the current hour, or left the location during the hour (Table 1). Action features included the participant's current activity that could be performed in those states (Table 2). Participants could stay for another hour, leave the location, and once they have left a location go to another location. The data contained a total of 149 days.

We modeled the state transition probabilities ( ) as a stochastic MDP to model the environment's influence on arrival time to a destination. The participants could stay or leave a place with 100% probability. Once the participants leave their current location, their arrival time at their destination depends on their desired arrival time and the environment (e.g., traffic, travel distance). To model the influence of the external variables, we empirically estimate the probability that participants have arrived at another place within an hour or not. The median number of states and actions per family were 14,113 and 85 respectively, for all combinations of possible features.

Aggressive Driving Behavior Data Set

Drivers that routinely engage in aggressive driving behavior present a hazard to other people in traffic [START_REF]Aggressive driving: Research update[END_REF]. To understand aggressive driving routines, it is important to explore the types of contexts aggressive drivers are likely to prefer (e.g., turn types, car speed, acceleration) and the driving actions they apply in those contexts (e.g., throttle and braking level, turning) (RQPat). Aggressive drivers might also be prone to dangerous driving behavior that does not occur frequently (e.g., rushing to clear intersections during rush hour [START_REF] Shinar | Aggressive driving: an observational study of driver, vehicle, and situational variables[END_REF]). Such behavior might manifest itself as variations from established routines (RQVar).

It is also important to compare the routines of aggressive drivers with non-aggressive drivers to understand how aggressive drivers can improve their routine (RQComp). To understand those differences, it is not enough to compare the contexts both groups of drivers find themselves in, but also the actions that drivers perform in those contexts. This is because both aggressive and non-aggressive drivers can attain similar driving contexts, but the quality of the execution of driving actions may differ. For example, both types of drivers might stop at a stop sign on time, but aggressive drivers might have to brake harder or make more other unsafe maneuvers than non-aggressive drivers.

This data set contains driving data from 22 licensed drivers (11 male and 11 female; ages between 21 and 34) from a mid-sized city in North America [START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF]. Participants were asked to drive their own cars on their usual daily driving routes over a period of 3 weeks. Their cars were instrumented with a sensing platform consisting of an Android-based smartphone, On-board Diagnostic tool (OBD2), and an inertial measurement unit (IMU) mounted to the steering wheel of the car. Ground truth about participants' driving styles (aggressive vs. non-aggressive) was established using their self-reported driving violations and responses to the driver behavior questionnaire [START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF]. The driving data collected in the study included: car location traces (latitude and longitude), speed, acceleration, engine RPM, throttle position, and steering wheel rotation. Sensor data was recorded every 500 milliseconds.

We use a subset of this data focused on intersections (where instances of aggressive driving are likely to occur [START_REF] Shinar | Aggressive driving: an observational study of driver, vehicle, and situational variables[END_REF]). We used location traces of the participants' driving routines to manually label intersections and the position of the vehicle in those intersections. One of the limitations of this data set is that there is no information about other vehicles and traffic signs and signals that represent the environment. We then split the intersection instances into sequences of sensor readings that start 2 seconds before the car enters the intersection, and end 2 seconds after the car exits the intersection. This resulted in a total of 49,690 intersections from a total of 542 hours of driving data from 1,017 trips.

To model states we combined the driver's goals (e.g., make a right turn), the environment (e.g., position in intersection), and the current state of the vehicle (e.g., current speed) into features of the states (Table 3). Actions in our model represent how the driver operates the vehicle by steering the wheel, and depressing the gas (throttle) and brake pedals. We aggregate the driver's actions between different stages of the intersection and represent the median throttle and braking level, and note any spikes in both throttle and braking. We consider the movement of the steering wheel to estimate whether the driver turned in one smooth action, or if the turn required one or more adjustments. Table 4 shows action features in our model. We identified 7,272 different states and 446 different actions in the data set.

VALIDATING THE MODELS OF ROUTINE BEHAVIOR

In this section, we evaluate the quality of the routines extracted using our model from the two data sets. First, we show that the routine actions we extract are predictive of the majority of behaviors in the data; i.e., that the algorithm is sufficiently predictive for modeling routines. Accuracy of this prediction task also quantifies the variability of the routines in the model, where high accuracy suggests low variability. It also shows that the extracted routines generalize to contexts and actions that we have not observed during model training. Second, we show that the routines extracted using our approach are meaningful. We show that the patterns extracted using our approach correspond to the actual routines and routine variations in our two example behavior logs (RQPat & RQVar). We also show that those routines and variations show the real differences between modeled populations (RQComp).

Quantifying Routineness of Human Behavior

Although we are not interested in the predictive power of the MaxCausalEnt IRL per se, we use the task of predicting the next action given a state to evaluate our model's ability to extract routine. Using 10-fold cross validation for each person in each dataset, we compare the performance of this algorithm for extracting routine behavior with a simple Zero-R algorithm, which always predicts the overall most frequent action, and a first-order Markov Model algorithm, which always predicts the most frequent action for each state. We chose these two baselines because they explicitly establish the frequency of actions in the training set. Matching or exceeding these baselines means that the algorithm has correctly identified frequent routine actions and that the predictive power of the algorithm is sufficiently high to model routines. We now show that the routine patterns extracted using our approach match the actual routines of people. To do this, we recruited researchers that work with machine learning and data mining in the domain of human behavior, and asked them to identify the routines and variations extracted using our approach. We then confirmed that those patterns matched the ground truth behaviors established in the previous work [START_REF] Davidoff | How routine learners can support family coordination[END_REF][START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF]. This allowed us to verify that the patterns extracted using our approach are meaningful and represent the actual routines.

Study Software

To make the routine behavior models created using our approach accessible to participants and allow them to investigate the extracted routine patterns, we developed a simple visualization tool. To maintain a level of familiarity, we base our visual encoding of routine behavior elements on a traditional visual representation of an MDP as a graph (Figure 1). Our MDP graph contains nodes representing states (as circles) and actions (as squares), directed edges from state nodes to action nodes (indicating possible actions people can perform in those states), and directed edges from actions to states (indicating state transitions for any given state and action combination).

To enable participants to see changes in features of states and actions, we encode state features and action features as a series of color-coded circular marks arranged in a spiral shape within the nodes. Each feature has a dedicated hue. Feature values that are present in the node are represented by a dark shade, and feature values not present in a light shade of that color. A dark boundary serves as a separator between features. More details, in text are always available simply by moving the cursor over a node (Figure 1.C).

To show frequent behaviors in the model, we visually represent the probability of different graph elements using line thickness. Thickness of the outside line of the state and action nodes encodes the frequency of that state in a behavior sequence ( ), where thicker lines indicate states that are likely to be part of a routine. Similarly, the thickness of the edges encodes the probability of that edge. Thickness of edges from states to actions is given by the probability distribution of actions given states ( ), and represents the influence of each state on the choice of actions. The thickness of edges from actions to states is given by the probability of transition ( ).

To layout the nodes, we sort the initial states from the demonstrated sequences by their frequency ( ) in descending order. We then use a version of the depth-first search algorithm, starting from the initial state nodes, that traverses nodes by first iterating over edges in order from highest to lowest probabilities and . State nodes are never duplicated (i.e., there is exactly one node in the layout for each state in the model), whereas action nodes are duplicated for each state.

When the participant selects a data set and population, the tool provides the initial layout of routines extracted using our algorithm. This shows the most important information about the extracted routines. However, to further analyze the routine behavior, the participant must be able to explore the details of routine variations filtered out as aggregate nodes (Figure 1.B.3). For example, the participant might want to find which of the parents' routines include locations where they pickup and drop off their children.

Aggregated items contain valuable information about potential routine variations. For example, a child might go to her grandparents or their friend's house on Wednesdays after school; two variations on the same routine that occur with similar probability. To show all possible variations in an aggregate, the participant can click on the aggregate to expand its content. To mark an aggregated node as a variation of interest, the researcher can pin that aggregated node by clicking on it, thus removing it from all of its aggregate parents. When the researcher holds Alt-key and clicks on an aggregated node, this pins all the nodes on the most likely sequence of states and actions, determined by the probabilities of edges between the two. The pinned sequence, starting from the clicked node to the sequence end node, represents a routine variation. Pinned nodes are identified by a gray glow effect. All nodes that are part of the extracted routine are automatically pinned, and all other nodes are unpinned. Clicking on a pinned node unpins it, which returns the node into the aggregate.

To determine whether or not to pin the node, the researcher can review the features of individual and aggregate nodes by hovering over them. In addition to showing the details of individual nodes in the details panel (Figure 1.C), hovering over nodes shows relationships between different elements of the routine. Hovering over any node highlights the most likely routine path from an initial state to the collecting node that contains the hovered over node (Figure 1.B.1). This makes it easier to understand the routine states and actions in the area of interest.

Participants

To verify the routines, we have recruited 8 researchers (5 male and 3 female) that have had experience with machine learning and data mining, or have worked in the domain of activity recognition and human routine modeling. The participants included Ph.D. students, Postdoctoral fellows, Research Scientists, and Professors working or visiting at our local University. All participants had experience with machine learning, 2 with data mining, 4 with activity recognition, and 1 worked specifically on modeling human routine behavior. The participants were compensated $25 for taking part in the study.

Method

When participants arrived at our lab, we briefed them on the purpose of the study and they signed a consent form. They then filled out a short questionnaire asking them about their occupation and experience with relevant research topics. We then demonstrated the visual tool to the participants and allowed them to practice using it for approximately 20 minutes. Participants had to complete two tasks: 1) identify daily routine for a randomly chosen person and weekday from the daily routines data set, and 2) identify the differences between routines of non-aggressive and aggressive drivers from the driving data set. The first task took approximately 20 minutes, and the second task about an hour. Total study time was approximately between one and a half and two hours.

For the two tasks, we asked participants to identify routines and differences between the routines without presenting them with the ground truth. We did this to avoid biasing the participants towards trying to match the routines presented in the tool with what we might have told them is the correct answer. We then compared their answers with the ground truth from the previous literature to verify that the tool extracted the right and meaningful routine. Because the main purpose of the study was to validate the extracted routines and not evaluate the usability of the tool, participants could ask clarifying questions about the tool and the user interface at any point during the tasks.

Ground Truth

In the first task, we used the family daily commute routines model to understand the patterns of pickups and drop-offs. We compare the findings of our participants with the ground truth and discussion provided by Davidoff et al. [START_REF] Davidoff | How routine learners can support family coordination[END_REF]. In this data set, the ground truth represents selfreported daily commute routines for all family members in all families that took part in the study. Family members reported their location and the time they usually arrive and leave that location. Davidoff et al. [START_REF] Davidoff | How routine learners can support family coordination[END_REF] then manually annotated and confirmed the routines in the raw sensor data.

In the second task, we used the driving routines of aggressive and non-aggressive drivers that Hong et al. [START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF] identified in their data set. Hong et al. [START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF] used their intuition and expert knowledge of driving behaviors to separately compare the distributions of each sensor stream in the raw data to gain insight about aggressive driving styles. They found that aggressive drivers drive faster than non-aggressive drivers, and that they experience higher acceleration than non-aggressive driver (i.e., they are more likely to press hard on the gas and brake pedals). Additionally, they found more variability in the behavior of aggressive drivers than non-aggressive drivers.

Results

Our results show that our approach extracts meaningful patterns of routine behavior. The participants were able to point out the patterns that form the high-level routines present in the ground truth for both tasks (RQPat). For the daily routines task, this means that they successfully listed the locations and times of the routines of people in the daily routine data set. However, two participants identified two separate patterns where the locations and times reported as part of the main daily routine did not correspond to the ground truth. After careful examination, we found that the participants wrongfully identified the actions as part of the routines because the events were part of infrequent routine variations that the people in the original study [START_REF] Davidoff | How routine learners can support family coordination[END_REF] did not report in the ground truth. The algorithm correctly assigned low probabilities to those actions, but the participants did not notice this. This is likely an issue with the visualization rather than the algorithm itself, and is something that can be addressed with more training with the tool.

In addition to simply pointing to patterns that represented correct routines, participants also generated some insights for themselves. For example, all six participants that were presented with a parent's daily routine that contained a child pickup or drop-off specifically pointed out this activity. Also, three participants, that had a case where the parent drops off the child as part of his or her routine, but does not also pick the child up, correctly explained that the other parent was likely responsible for the pickup, without seeing the other parent's routines.

In the driving data set, participants pointed to the patterns that form the main routines (RQPat) and variations in driving behavior (RQVar) of both aggressive and nonaggressive drivers. All participants pointed to patterns that show that aggressive drivers are more likely to drive faster through intersections than non-aggressive drivers. Five participants showed the patterns of routine variations where aggressive drivers are likely to increase their throttle just before entering and leaving intersections. Participants pointed those out as the main differences between the two populations (RQCom). Two participants also pointed to the probabilities of routine variation patterns extracted using our approach that suggest that aggressive drivers are less consistent in their behavior than non-aggressive drivers. Participants likely drew their conclusions from the model, but might also have a preconceived notion that acceleration and speed are correlated with aggressive driving. However, even if our participants had preconceived notions, they could verify and document them using our model.

Although evaluation of the visualization tool was not our main goal, 2 participants mentioned that such a tool would help them explore and understand human routines. One participant, who studies routine behaviors, pointed out that the organization of routine patterns and variations helped him clarify his understanding of what constitutes routines and how they manifest themselves in people's activity.

DISCUSSION

Through our evaluation, we showed that our models trained using MaxCausalEnt algorithm [START_REF] Ziebart | Modeling interaction via the principle of maximum causal entropy[END_REF] can extract patterns of routine behavior from demonstrated behavior logs. This is a novel application of an algorithm that was designed to predict human behavior. In our evaluation of the algorithm we found that its ability to predict routines from the two example data sets was sufficient for modeling routines. The ability of MaxCausalEnt algorithm to generalize from small sample sizes enabled it to beat the baseline in the daily routine data set. The performance of the algorithm was comparable with the first-order Markov Model in the aggressive driving data set. This is likely because the training data happened to match the testing data well. However, this is not safe to assume in general case, and MaxCausalEnt's decision-theoretic guarantee that it will not overfit the observed data make it a better choice for modeling routines than the first-order Markov Model.

Although not our main contribution, our visualization tool also helped us validate the ability of our approach to extract meaningful routines. The participants were able to explore context, actions, and the relationships between the two, to correctly identify the patterns of routines (RQPat) and their variations (RQVar) reported in the previous work. They pointed to these relationships to establish the differences between routines of the two driver populations (RQComp). Although we carefully designed our tool, our goal was not to formally evaluate its usability. We did not notice any usability issues that prevented participants from learning the elements of the model. We found that the participants knew how to progress towards understanding the routines.

Our results imply that researchers can use the patterns extracted using our approach to more quickly identify major aspects of routines by visually inspecting them, even after only short amount of training, compared to previous work. For example, Davidoff et al. [START_REF] Davidoff | How routine learners can support family coordination[END_REF] performed tedious manual labeling of routine and routine variation patterns in the raw data based on feedback from the participants before presenting the patterns on a timeline. Hong et al. [START_REF] Hong | A smartphone-based sensing platform to model aggressive driving behaviors[END_REF] used their intuition and expert knowledge of driving behaviors to separately compare the distributions of each sensor stream in the raw data to gain insight about aggressive driving styles. Our participants had to only explore the patterns extracted using our approach.

The knowledge that the researchers gain about routine behaviors through exploring our models can inform the design of interventions that help people improve their routines. For example, the knowledge that aggressive drivers are likely to use higher throttles can inform the design of in-car systems that monitor the throttle and make the driver more aware of this aggressive behavior through subtle ambient notifications. Another advantage of our approach is the underlying MDP-based model, which can be used to power smart agents that automatically classify current behaviors and prescribe new actions that improve existing routines.

CONCLUSION AND FUTURE WORK

We presented a novel approach for modeling human routine behavior from behavior logs that explicitly models one of the most important aspects of routines: the causal relationship between the contexts and actions the people perform in those contexts. We demonstrated that our approach can be used to extract meaningful routine patterns from two different types of human behaviors. However, future work should explore how our approach can be used to model routines in other domains (e.g., health, accessibility, software user interfaces).

We showed that routine models extracted using our approach can help researchers to identify routines and routine variations without having to manually search for those patterns in raw data. This is a step towards models that enable researchers to form scientific insight and generate knowledge about human behavior that will help them develop new theories about human routines. Having such knowledge can also inform the design of novel technologies that help people improve the quality of their routine behavior.

Our long-term goal is to help researchers in the process of generating knowledge that will support the design of technologies for routine behavior change. In this work, we presented a simple visual representation of the model extracted using our approach. However, in future work, we intend to develop a visual analytics [START_REF] Daniel | Challenges in Visual Data Analysis[END_REF] tool to offer a framework for visual data mining and understanding of human routine behaviors. We also intend to show how researchers and designers can use the insights they generate from our models to design technologies that help people improve their routines.

Machine understanding of human routine behaviors is also necessary to design and implement effective technologies that help people improve their routines. Our approach encodes patterns of routine behavior in a way that allows systems, such as smart agents, to classify, predict, and reason about human actions under the inherent uncertainty present in human behavior. For example, future smart agents can use our model to detect aggressive driving behavior and provide feedback to the driver based on the models of non-aggressive drivers. Such technologies can have a positive effect on society by making people healthier, safer, and more efficient in their routine tasks.

Figure 1 .

 1 Figure 1. Study software user interface showing the main routine and one likely variation of non-aggressive drivers extracted using our approach: A) overview panel, B) the main display area containing subgraphs representing automatically extracted routine sequences of states (circles) and actions (squares), B.1) the user is hovering over a node to highlight extracted routine (nodes highlighted in purple and dark gray edges), B.2) an action that starts a variation from the main routine B.3) aggregate items representing possible extracted variations, and C) details panel showing information about visual elements on demand.

Table 1 . State features capturing the different contexts of a daily commute.

 1

	Feature	Description
	Day	Day of week {M, T, W, Th, F, Sa, Su}
	Time	Time of day in increments of 1 hour {0-23}
	Location	Current location
	Activity	Activity in the past hour
		{STAYED AT, ARRIVED AT, TRAVELING FROM}

Table 2 . Action features representing actions that people can perform when at a location.

 2

	Feature	Description
	Activity	Activity people can perform in current context
		{STAY AT, TRAVEL TO}
	Location	The current location to stay at or next location to go to

Table 3 . State features capturing the different contexts the driver can be in.

 3

	Feature	Description
		Goals
	Maneuver	The type of maneuver at the intersection
		{STRAIGHT, RIGHT TURN, LEFT TURN, U-TURN}
		Environment
	Position	Current position of the car in the intersection
		{APPROACHING, ENTERING, EXITING, AFTER}
	Rush hour	Whether the trip is during rush hour or not
		{TRUE, FALSE}
		Vehicle
	Speed	Current speed of the vehicle (5-bin discretized)
	Throttle	Current throttle position (5-bin discretized)
	Acceleration	Current positive/negative acceleration (9-bin discretized)
	Wheel Position	Current steering wheel position
		{STRAIGHT, TURNING, RETURNING}
	Turn	Current turn vehicle is involved in
		{STRAIGHT, SMOOTH, ADJUSTED}

Table 4 . Action features representing actions that drivers can perform between stages of the intersection.

 4

	Feature	Description
	Pedal	Median throttle (gas and brake pedal) position
		(10-bin discretized)
	Throttle Spike	Sudden increases in throttle
		{NONE, SUDDEN, INTERMITTENT}
	Brake Spike	Sudden braking
		{NONE, SUDDEN, INTERMITTENT}
	Turn style	Type of turn driver performed in intersection
		{STRAIGHT, SMOOTH, ADJUSTED}

ACKNOWLEDGMENTS

This work was funded by NSERC (PGSD3-438429-2013) and NSF (CCF-1029549, IIS-1217929). The authors would like to thank Brian Ziebart for his valuable input regarding the MaxCasualEnt algorithm, Scott Davidoff and Jin-Hyuk Hong for their help in obtaining and understanding the two datasets used in our work, and Julian Ramos and Christine Bauer for discussions about human routine behavior.


Pierre-Alexis Herrault

David Sheeren

Mathieu Fauvel

Claude Monteil

Martin Paegelow

A COMPARATIVE STUDY OF GEOMETRIC TRANSFORMATION MODELS FOR THE HISTORICAL 'MAP OF FRANCE' REGISTRATION

Keywords: Historical map registration, geometric correction, kernel regression

It is widely recognized that present-day biological diversity may reflect past land use and landscape features. Understanding the historical background is essential to explain the ecosystem functioning of landscapes today. Numerous historical spatial data have been used to reconstruct these pathways but resulted in problems of sources, formats and supports. One major problem is how to incorporate historical data in current coordinate systems to enable them to be compared with each other and with current sources. Previous works on the registration of historical maps already demonstrated the superiority of local methods (e.g. a Delaunay-based method) over global methods (like polynomial mapping models) because of the existence of local geometric distortions. The same studies also highlighted the importance of selecting ground control points with a homogeneous spatial distribution to improve the accuracy of registration. Furthermore, while kernel-based methods have already proved their efficiency for many other applications, they have rarely been used for map registration even though they provide an interesting alternative to conventional methods. In this paper, we present a comparative study of various geometric transformation methods to register an excerpt of the historical 'Map of France' (an Ordnance Survey map) dating from the 19th century. We compare the performance of several global and local methods with kernel-based methods (Gaussian and polynomials) and analyze the impact of the number of control points, their nature and their spatial distribution on the quality of registration. A protocol was developed to apply the transformation models in various situations and to identify the best strategy for the georeferencing of historical maps.

INTRODUCTION

It is widely recognized that present day biodiversity may reflect past land use because of a possible delay between changes in a landscape and changes in the associated plant or animal population. Knowing the historical background of landscapes is therefore a precondition to understanding ecosystems functioning [START_REF] Jackson | Vegetation, environment, and time: The origination and termination of ecosystems[END_REF]. Historical ecology concerns all the elements in a landscape. For example, using archival and archeological data (Hewlett, 1977), reconstructed river valleys in Great Britain from the year 1000. Later, other authors used similar approaches to explain current forest biodiversity with respect to forest history [START_REF] Hermy | An indirect gradient analysis of the ecological relationships between ancient and recent riverine woodlands to the south of Bruges[END_REF].

Different problems arise when users want to use historical spatial data, for instance maps drawn in the 19th century. Firstly, the variety of sources with different formats and supports increase the heterogeneity of spatial data. Secondly, because of the technical and financial means at that time, the quality of historical data may differ from one source to another and may lead to inconsistencies and biases in the assessment of landscape changes (Pontius and Lippit, 2006). Finally, these data are generally not georeferenced. In order to incorporate them in current coordinate systems and to compare them with present data, a registration step (or georeferencing) is required. Georeferencing is a crucial step in the use of historical data. Several studies have already focused on this issue (e.g. Baiocchi and Lelo, 2005;[START_REF] Bitelli | Ancient map comparisons and georeferencing techniques : a case study from the Po River Delta (Italy)[END_REF].

Generally speaking, two main approaches are used for map (or image) registration: global or local [START_REF] Zitova | Image registration methods: a survey[END_REF].The global approach is appropriate when map distortions are homogeneous. With this method, only one transformation model is computed and applied using all available ground control points (GCPs). However, a local approach is more efficient if local distortions have to be taken into account. With this method, the map or image is partitioned into a set of regions and one transformation model is computed for each region individually.

The local approach is known to be more efficient to register historical maps [START_REF] Fuse | Vizualizing the landscape of old time Tokyo (Edo City), Commission V, WG V/6. Goshtasby A[END_REF]. The low quality of the survey instruments, the lack of cartographical specifications and the diversity of the surveying engineers and campaigns led to numerous local planimetric errors that cannot be handled correctly by global transformation models. However, other registration techniques could be used to improve the accuracy of the geometric transformation models. In particular, nonparametric methods like kernel smoothing regression are a possible alternative to conventional approaches because they are locally sensitive but provide one global model. These methods have already proven to be very efficient for many applications including image processing and pattern recognition (Shawe-Taylor and Cristanini, 2004;[START_REF] Takeda | Kernel regression for image processing and reconstruction[END_REF].

Like the number of points, the spatial distribution of the GCPs is also important in the map registration process since it may affect the accuracy of the transformation models [START_REF] Mather | Map-image registration accuracy using least-squares polynomials[END_REF][START_REF] Wang | Effect of the sampling design of groundcontrolpoints on the geometriccorrection of remotely sensed imagery[END_REF]. Regular spatial distribution of GCPs provides a bestfitting process but requires many points for efficient registration [START_REF] Boutoura | Some fundamentals for the study of the geometry of early maps by comparative methods[END_REF]. In the case of old topographic maps, the GCP problem is often difficult to overcome due to the possible change in the landscape objects over the past centuries and in the characteristics of the GCPs themselves (reliability and accuracy, kind of object, graphic representation). In contrast to aerial and satellite images, the effect of the sampling design of the GCPs, their number and their category on the geometric correction of old maps have been less studied.

The aims of this paper are to test and assess the potential of the kernel regression method to register an historical map from the 19th century. We compare the method with other global and local conventional approaches. Second, we independently investigate the effect of three parameters used to select GCPs on the accuracy of map registration:

1. The spatial distribution of the GCPs (regular or random); 2. The number of GCPs; 3. Whether the GCPs are reliable or unreliable.

The paper is organized as follows. In section 2, we describe the data and the experimental protocol with the kernel regression and other registration methods. In section 3 we present and discuss the results of the experiments. Finally, some conclusions are drawn along with some perspectives.

MATERIAL AND METHODS

Data

In this comparative study, we used an excerpt of the historical 'Map of France' produced from 1825 to 1866 at 1:40000 in the projection of Bonne. This excerpt covers an area of 20 km² in southwestern France surrounded by contrasting landscapes: valleys bottoms, mountains and a plateau. This historical Map of France is known for its planimetric accuracy, considerably better than the previous Cassini map published during the reign of Louis XV. Moreover, the Map of France is in color, which makes it easier for users to capture the cartographic elements. Last, the Map of France was produced in an interesting period for the study of the forest evolution. The extent of French forests was lowest around 1830 and then expanded continuously until today [START_REF] Mather | The course and drivers of the forest transition: The case of France[END_REF]. Thus, if an existing forest appears on the old map, it might be unchanged for several decades. The GCPs were selected from an existing database of aerial photographs: the BDOrtho® produced by the French National Mapping Agency (IGN). The image has a spatial resolution of 0.50 m.

Methodology

An experimental protocol was defined (Fig. 2) to compare the registration methods. The protocol comprises four main steps: (1) acquisition of a high number of two categories of control points in the old and reference maps according to two spatial distributions (random and regular), (2) creation of GCP subsets of various sizes, (3) computation of several transformation models, (4) validation and comparison of the results. In the first step (1), we identified landmarks and crossroads to assess the impact of the categories of GCP on the accuracy of the geometric transformations. We assumed that the landmarks are more accurate than the crossroads. Indeed, these points (chapels, churches, mills, etc.) were used as reference points when topographical surveys were being conducted before the map was drawn. For this reason, their absolute planimetric accuracy might be higher than other points on the map, particularly the crossroads.

A total of 75 randomly distributed landmarks were found on the old map. For the crossroads, we selected 150 GCPs according to two spatial sampling processes: random and regular. The regular distribution was obtained from a regular grid covering practically the entire map and divided into 12 x 20 squares. For the landmarks; we only used random distribution because the small number of points prevented regular sampling.

In the next step (2), we created GCP subsets with the aim of assessing the impact of the number of points on the accuracy of the transformation models. The number of points in these sets ranged from 25 to 150 for the crossroads and from 25 to 75 for the landmarks.

In step (3), we compared the performance of several geometric transformation models for each configuration of the GCPs. Global and local methods were tested and compared to the alternative kernel regression approach. The methods tested are listed in

• Conventional global and local registration methods

As mentioned above, previous works on the registration of historical maps have showed the advantage of local methods (like the triangle-based method) over global methods (polynomial mapping functions) because of the existence of local geometric distortions (White and Griffin, 1985;[START_REF] Fuse | A study on geometric correction of historical maps[END_REF].

A global transformation simply means that all the control points are used to derive a single mathematical model of the geometric distortions. First order, quadratic or cubic models are the most widely used to register an image or a map. The general expression of the 2D polynomial model is defined as follows:

1 0 0 N N i j ij i j x a u v ¦¦ 1 0 0 N N i j ij i j y b u v

¦¦

where a ij and b ij are the coefficients of the equations ݔ = ,ݑ(݂ )ݒ and ݕ = ,ݑ(݃ ,)ݒ the reference coordinates being ,ݑ( .)ݒ The degree of the polynomial model is defined by ܰ. For ܰ = 2, the polynomial function can be expressed as follows: A minimum of three points (6 parameters) or six points (12 parameters) are needed to estimate the coefficients of the first order model (linear transformation) and second order model (nonlinear transformation) respectively. However, the number of equations is usually higher than the number of unknowns (i.e. there are more GCPs than the minimum number required). The coefficients of the mapping functions are then computed using least square fit, which provides residuals for each point.

In the case of the local transformation approach, after dividing the map into a set of regions, several local mapping functions are computed (Goshtasby, 1986[START_REF] Goshtasby | Piecewise cubic mapping functions for image registration[END_REF]. In general, the regions are created from a Delaunay triangulation and one first order polynomial function is computed per triangle. Each triangle is based on three GCPs and define a mapping function that is only valid within this region. According to (Balletti, 2006), the smaller the triangular areas of the mesh, the better the results. Thus, for this approach, it is important to collect a large number of control points to reduce the extent of the triangle.

• Kernel-based method

In this paper, we compare the performances of several global and local methods with kernel-based methods (Table no. 1). These models belong to the category of global methods but are able to take local distortions into account by combining local models computed in a bandwidth (locally weighted regression smoother, Fig. 3). (from [START_REF] Takeda | Kernel regression for image processing and reconstruction[END_REF]))

The previous conventional methods were compared with the kernel regression (Shawe- Taylor and Cristanini, 2004). Kernel regression is a non-parametric statistical technique that belongs to the category of global methods but which is locally sensitive since it enables local distortions to be taken into account by combining local models computed in a given bandwidth (locally weighted regression smoother). This method is efficient for the estimation of nonlinear relationships between variables with minimal assumptions about the model. The general expression of the model is given by the decision function (given here for the 1-D case for the sake of simplicity):

)ݖ(݂ = ߙ ݔ(݇ , )ݖ + ܾ ୀଵ
where ݇(. . , . . ) is the kernel, ߙ the amplitude coefficients, ܾ the bias, ݔ the learning samples and z, the predictive sample. The ߙ and ܾ are the parameters to optimize. They are determined by minimizing the mean square error (MSE), i.e.:

݉݅݊ ఈ , (ݕ െ ݔ(݂ )) ଶ + צߣ צ݂צ צ ଶ ୀଵ
where the ߣ parameter is a regularization parameter that is used to prevent from overfitting.

After straightforward computations, the solution is given by (in matrix form):

ߙ ො ܾ ൨ = ܭ(ݕ + ܫߣ ାଵ ) ିଵ
where )ݕ( = ݕ , ܭ( ) , = ݔ(݇ , ݔ ܭ(,) ) ାଵ, = ܭ( ) ,ାଵ = 1, ܭ( ) ାଵ,ାଵ = 0 and ܫ ାଵ is the identity matrix of size n+1.

The kernel is a positive semi-definite function. It is used to map the sample implicitly in a feature space where the relation between X and Y is easier to estimate.

In this paper, two kernels were investigated:

-Polynomial kernel: ܭ = ݔۃ൫ , ݔ ۄ + 1൯ -Gaussian kernel: ܭ ఙ ,ݔ( )ݕ = exp ቀെ ԡ௫ି௬ԡ మ ଶఙ మ ቁ
In the final step of the protocol (step 4), the results were validated and compared to identify the best registration approach. We used cross validation to assess the accuracy of each transformation process and in particular, the leave-one-out method. Cross validation is commonly used to estimate the accuracy of a predictive model. It involves dividing the observations into two subsets: one subset is used to train the predictive model and the other to validate it. This process is repeated several times to increase the reliability of the validation.

In the leave-one-out cross-validation, only a single observation is selected for validation and all the remaining observations are used as training data. This process is repeated such that each observation in the sample is used once as the validation data. The validation results are then averaged to compute an average error. In our case, the error index used for assessing the quality of the registration was the root mean square error (RMSE). Thus, after the cross-validation process, we obtained an average value of the RMSE (called RMSEm).

RESULTS AND DISCUSSION

We successively analyzed each parameter of the registration process in order to identify the best strategy for georeferencing the historical 'Map of France'.

• Impact of the number of GCPs

This parameter was only assessed for the crossroads, of which there were many examples. Table no. 2 summarizes the results of different geometric transformation models using the five GCP subsets composed of 25, 50, 75, 100 and 150 points. The influence of the number of points was analyzed for the two spatial distributions available: random and regular.

Overall, the results showed that the number of errors decreased with an increase in the number of points regardless of the spatial distribution. For instance, we observed that, for the affine transformation, the RMSEm decreased by approximately 27 m from the GCP subset of 25 points to the subset of 150 points. The same was true for the kernel transformation. Using the Gaussian kernel, the average error decreased by 26 m (in the case of regular distribution) and by 37 m (in the case of random distribution) using a number of GCPs that ranged from 25 to 150. The Delaunay-based transformation was more sensitive to the number of points. The RMSEm followed the same negative trend but the initial number of errors was much higher than with the other transformation models. When the number of the GCPs is small, the size of the triangles is very large. Consequently, the local first order polynomial models are less accurate.

Table no. 2. RMSE(m) obtained by the different transformation models tested for 25 to 150

GCPs and for a regular or a random distribution of these GCPs.

• Impact of the spatial distribution of GCPs

Concerning spatial distribution, Table no. 1 shows that the regular distribution is generally more efficient than the random one up to a set comprising 75 GCPs. While the error for the affine transformation varied between 131 and 88 ms for a random sampling, the average error ranged from 106 to 6 meters when the spatial distribution was regular. This is also true for the first and second order polynomial kernels. However, with over 75 GCPs, the trend was reversed in favor of random distribution (Fig. 4). The most accurate models are computed from a random sampling when large GCP datasets are used. In the case of the Gaussian kernel, the random distribution is always better. • Impact of the geometric transformation models

The performances of the geometric transformation models were only compared for the random distribution (Fig. 5).

Generally speaking, our results show that the differences in accuracy between the models are greater with small sets of control points than with large sets. When only 25 points were used, the best model was given by the Gaussian kernel (RMSEm equals 75 m). The least efficient model was the Delaunay-based transformation (RMSEm equals 195 m). As mentioned above, this local approach appears to be very sensitive to the number of points.

The results of the Gaussian kernel model were always better than the results of the other methods (except for the set of 75 GCPs). The Gaussian kernel model appears to be very efficient when the number of GCP is high (RMSEm equals 42 m for 150 GCPs). The third order polynomial kernel also gave good results with large datasets with an average error equal to 52 meters for 150 GCPs, up to 10 meters more than the Gaussian kernel transformation. Kernel-based methods are generally more efficient than conventional global models for equivalent transformation (see Fig. 5). This was the case for almost all the GCP sets (except for the set of 100 GCPs). • Impact of the GCP categories Finally, we assumed that landmarks (or remarkable points) are more accurate and provide better results.

Globally, this assumption was confirmed. In the majority of the cases, the RMSEm is lower when landmarks are used in geometric transformation models.

With the affine transformation the difference in errors between landmarks and crossroads ranged from 33 to 16 meters. The same relation was found for second order polynomial transformation and for kernel based transformations in general. For example, while Gaussian kernel transformation presents an average error of 64.12 meters with crossing roads (set of 75 GCPs), the mean of positional error was about 59 meters with landmarks for the same number of GCPs.

The trend was different with the second order polynomial and with the Delaunaybased transformation.

With the second order polynomial model, there was a slight difference in the number of errors using landmarks and crossroads, but overall, crossroads produced better results. The difference ranged from 2 to 4 meters. This opposing trend can be explained by the random sampling process. This step was performed using a GCP with 200 crossroads but only 75 landmarks. Furthermore, the standard error of the order three polynomial model using crossroads was much higher than that using landmarks so the variability of the results was not the same, hence explaining these contradictions. As regards the Delaunay-based transformation, as we mentioned above, this model is very sensitive to the number of points. This time, the random sampling process using the crossroads category led to better spatial distribution. The probable explanation is that the size of the triangle was more regular in a few samplings, thus avoiding abrupt bias faults in some parts of the map.

Fig. 6 compares the transformations between the crossroads and the landmarks for 75 GCPs. Except for the Delaunay transformation and the third order polynomial, the landmark sets gave better results.

• Comparative analysis

The comparison of different geometric transformation models and the evaluation of some parameters in the selection of GCPs led to various teachings in order to find the best strategy.

First, the higher the number of GCPs, the lower the average number of errors. This is to be expected when predictive models are used for image registration: the higher the number of points used, the more local distortions are taken into account, and the more efficient the transformation will be. Moreover, the larger the sample, the lower the approximation so the model predictions are better and provide a good estimation of the transformation. Concerning this parameter, Delaunay triangulation appears to be a special case. Indeed, an increase in the number of GCPs leads to a spectacular decrease in the RMSE(m). This can be explained by the reduction in the size of the triangle that consequently takes local distortions better into account. With the other transformation models, the decrease in the number of errors between each GCP set was lower because the map registration is homogeneous and accomplished using all the GCPs.

The evaluation of spatial distribution revealed two major trends. When the number of GCPs is small, regular distribution is required. The explanation may be that regular distribution ensures global coverage of the region of interest. Even if the number is small, each part of the map is taken into account by the models, thus avoiding abrupt bias fault. A random distribution could lead to parts in the map being neglected and for that reason, not provide an accuracy transformation in all cases.

However, when the number of GCPs is high, random distribution is more efficient. Planimetric accuracy is not homogeneous in historical maps. Thus, local distortions can be solved by random distribution because the GCPs are not spatially constrained. Based on random sampling, random distribution appears to be more efficient for the registration of historical maps.

Next, we compared some global and local models and tested the potential of kernel based methods.

Overall, kernel based methods proved to be more efficient than conventional methods. When we focused on transformations with GCP sets numbering 150, the results of polynomial kernels were better than those of conventional polynomial methods, except for second order polynomials, whose results were very close.

The results of the Gaussian kernel method were the best. With 150 crossroads and random distribution, the average positional error was about 42 meters. This error is clearly lower than other transformations meaning that this method is a very accurate registration model and a very good alternative to methods currently employed in georeferencing of historical maps.

Finally, we analyzed categories of points. Distinguishing between landmarks and crossroads, we assumed that some points on the early map are reliable and others unreliable.

The results indicated that models using landmarks provided better transformations than the others. This result confirms that survey points are recommended in the creation of a GCP set. Nevertheless, their limited number (only 77 in this excerpt) meant we could not only use this category. Combining landmarks with crossroads was required to create GCP sets and to improve the transformations.

CONCLUSION

From this study, it appears that georeferencing of old maps depends not only on the transformation model used but also on the spatial distribution, the number, and the category of ground control points.

Kernel based methods proved to be more efficient than conventional methods. The with currently used methods showed that kernel models enable early maps to be georeferenced most accurately. In particularly, the Gaussian kernel provided the best transformation.

Analysis of the spatial distribution of GCP, their number and their category, highlighted an optimal configuration for georeferencing. Indeed, random distribution is required to better account for local distortions in the historical map. In addition, the models are more efficient with a large number of points for two major reasons: the region of interest is more comprehensively covered and error estimation by the prediction models is better. Finally, the analysis of the categories of GCPs showed that landmarks are more accurate than crossroads, but that their number is not sufficient for optimal transformation. This research could be extended to cover the following points:

First, the evaluation of Delaunay transformation showed that it was very sensitive to the number of GCPs. With 150 GCPs, it was less efficient than kernel based methods. However, with a large number of GCPs, it was more accurate than the other models tested . Combining landmarks and crossroads could thus greatly improve the results. Also, other categories of points could be analyzed to find the best combination of GCPs. Finally, we believe it is important to work on visualization of the error of transformation. This information could advance our understanding of local distortions to improve the local calibration of transformation models.

Fig. 1 .

 1 Fig. 1. Excerpt of the Historical 'Map of France' used in this study with examples of forest features.

Fig. 2 .

 2 Fig. 2. Experimental protocol applied to evaluate the effects of the nature, the number and the spatial distribution of GCPs on the registration accuracy.

  v b uv b u b v

Fig. 3 .

 3 Fig. 3. Local models successive (red lines) and the global kernel regression (blue line).(from[START_REF] Takeda | Kernel regression for image processing and reconstruction[END_REF]

Fig. 4 .

 4 Fig. 4. Comparison of the performance of registration methods for 150 GCPs and for a random or a regular distribution of these GCPs.

Fig. 5 .

 5 Fig. 5. Comparison of the performance of registration methods for 25 to 150 GCPs and for only one random distribution of these GCPs -Focus on the comparison between conventional and kernel registration methods.

Fig. 6 .

 6 Fig. 6. Comparison of the performance of registration methods for 75 GCPs according the nature of these GCPs.

Table no . 1.

 no

Table no . 1. The geometric transformation models compared. Global Method locally non sensitive Global Method locally sensitive Local Method

 no

	Affine	Linear Kernel	Delaunay
	Second order polynomial	Second order polynomial kernel
	Third order polynomial	Third polynomial kernel
		Gaussian kernel

Table no . 3. RMSE(m) obtained by the different transformation models tested for 25 to 75 GCPs with remarkable points (landmarks) or crossroads.

 no

Acknowledgments:

This research was supported by the French National Research Agency (ANR JCJC MODE-RESPYR 2010 1804 01 -01). P.-A. Herrault is also funded through a PRES Toulouse University and Region Midi-Pyrenees grant.

R E F E R E N C E S

Baiocchi V., Lelo K., (2005). Géoréférencement des plans historiques (du XVIIIe au XIXe) de la ville de Rome, et leur comparaison avec des cartes actuelles, Géomatique Expert 45, p. 42-47 Balletti C., (2006). Georeference in the analysis of the geometry content of early maps", E-Perimetron, 1-1: 32-42