C B S Tor

Lim

Z Sun

X P Tan
email: xptan@ntu.edu.sg

M Descoins

D Mangelinck

S B Tor

C S Lim

Revealing hot tearing mechanism for an additively manufactured high-entropy alloy via selective laser melting

Keywords: additive manufacturing, high-entropy alloy, laser deposition, hot tearing

Revealing hot tearing mechanism for an additively manufactured high-entropy alloy via selective laser melting

solidification cracking

Belonging to the category of powder-bed fusion techniques for metal additive manufacturing (AM), Selective laser melting (SLM) has the unique advantages of producing parts with refined surface finish, having a wide material library and fabricating complex geometries [1]. Besides these benefits, it is also proven recently to have the potential of manufacturing components with superior mechanical properties as compared to its conventional counterparts [2][3][4]. However, despite the continuous efforts worldwide, feedstock used in SLM are predominantly comprised commercial engineering materials [1].

Material development catered to metal AM is still limited which has inevitably hindered its industrial applications. Therefore, it is imperative to deploy new metals and alloys for SLM process. High-entropy alloy (HEA) is a relatively new class of alloys which normally has multiple principal elements, usually 5 or more, in equimolar or near-equimolar ratios [5,6]. Some HEAs have shown exceptional properties for cryogenic toughness [7] and strengthductility combination [8]. CoCrFeNi-based alloys have been one of the most widely researched HEA series [9][10][11][12]. Hence, it is of great interest for both academy and industry to achieve high-quality printing of HEA materials via SLM. Some CoCrFeNi based HEAs have been printed using SLM process [13][14][15][16].

However, it was found that under normal processing conditions, SLM-built CoCrFeNi exhibits poor tensile properties (e.g. ultimate tensile strength (UTS) is ~480 MPa and elongation is only ~ 8%) [13]. Carbon additions improved its UTS to ~800 MPa but the elongation was still below 15 %, as compared to the cast sample with an elongation of ~80% [9,14]. For SLM-processed CoCrFeNiMn, hot isostatic pressing (HIP) was adopted to improve the densification [16]. All these suggest that HEAs might be susceptible to hot tearing or other build defects which was not paid enough attention by the previous studies.

Lately, hot tearing phenomenon has increasingly been observed in metal AM. Though, most of the cases are confined to superalloys [17][18][19]. It is primarily due to the segregation induced precipitation and/or liquid film formation. Measures such as composition modifications and nano-particle additions were conducted to minimize the cracking [17,20]. However, it is still unclear (i) to what extent hot tearing happens in SLM-built HEA, and (ii) what causes it to happen since HEA usually has single-phase microstructure. Both above-mentioned questions are critical towards improving the printability of HEAs in metal AM [21].

Gas atomized pre-alloyed equimolar CoCrFeNi powder (~20 to 63 µm) was used on SLM 250 HL system (SLM Solutions, Germany). To investigate the three-dimensional (3D) printability of the material, a wide range of processing parameters were tested. Firstly, by keeping the laser power constant, energy density was varied from 50 to 180 J/mm 3 . 60 samples with various processing parameters were built. None of the samples has a relative build density higher than 99%. The optimal process parameter set holds a laser power of 150 W, hatch spacing of 100 μm, scanning speed of 270 mm/s and layer thickness of 50 μm. It generated a relative density of 98.7%. Secondly, two different scanning strategies, namely chessboard and stripe/bi-directional (illustrated in Fig. 1), with varying scanning lengths were deployed. For chessboard scanning strategy, laser scanning takes place randomly within the chessboard to minimize the residual stress. The representative SLM-built CoCrFeNi "C" sample (by chessboard scan) and "S" sample (by stripe scan) each has a strip size of 2.5 mm and 2 mm respectively.

Standard metallographic specimen preparation procedures for SLM-built stainless steel 316L (SS316L) were used on the HEA samples [22]. Optical microscopy (OM) and scanning electron microscopy (SEM) were carried out on Zeiss Axioskop 2 Mat and JEOL 7600F, respectively. Electron backscatter diffraction (EBSD) was performed on the AZtecHKL Oxford system with a step size of 3 μm. The misorientation analysis was conducted via the MTEX software [START_REF] Bachmann | Texture analysis with MTEX-free and open source software toolbox[END_REF]. Tensile coupons were machined via electron discharge machining (EDM) to a gauge length of 10 mm and a cross section of 3 x 2 mm 2 .

The tensile testing was conducted on an Instron static tester (series 5569) with a strain rate of 1 x 10 4 s -1 . Contact based extensometer was used for the test. At least 3 sets of tests were conducted for each scanning strategy. Micro computed tomography (μ-CT) images was acquired using the Skyscan 1173 system (Bruker, Kontich, Belgium). A high voltage of 130 kV was used to achieve a resolution of ~6.4 μm. Atom probe tomography (APT) analysis was performed by LEAP 3000X HR (CAMECA, Gennevilliers, France). Focused ion beam (FIB) on a FEI Helios dual-beam was used for sample lift-out. The micro-tips with an end radius of ~50 nm were prepared using the annular milling method. Pandat software package and existing solidification data were used to generate the solidification paths of stainless steel 316L (SS316L), CoCrFeNi and Al7075 [20,[24][25].

Fig. 1(a) and 1(b) display the cross-section images of the C and S samples along the build direction. The presence of cracks can be clearly noted. To reveal the detailed grain morphology, large-area EBSD scans with respect to the build direction were carried out. It is noted that the grain size of SLM-built CoCrFeNi samples is extraordinarily large, and columnar grains can epitaxially grow up to ~3 mm in length and ~200 μm in width. In comparison, SLM process typically generates a grain length of less than 0.1 mm for AlSi10Mg [26], and ~0.2 mm for SS316L, Inconel 625 and CoCrMo [4,27,28]. One close match is pure molybdenum with a grain length of ~2 mm, which also suffered from extensive intergranular cracking [29]. As known, lacking alloying content for grain nucleation will cause formation of big-sized grains in metals and alloys [30]. Thus, both CoCrFeNi and pure molybdenum might undergo limited constitutional undercooling during rapid solidification of AM process which leads to large columnar grains. Moreover, we found that epitaxial grain growth was broken down in the periphery of the lack-of-fusion cracks. However, small spherical pores circled in red have little influence on the grain growth. As observed by Leung et al. that spherical pores were mainly formed at the last stage of solidification due to the reduction of gas solubility as molten metals cooled down [31]. Therefore, epitaxial columnar grain growth is unlikely to be affected. Both samples exhibit a strong <011> crystallographic texture. By comparing the OM and EBSD images, ~70 cracks are randomly selected for misorientation computation. It is found out that all cracks are intergranular in nature and majority of the adjacent grains' misorientations lies in between 40° to 50° (Supplementary Fig. 1).

Tensile testing results show that the yield strengths (Ys) are 556.7 ± 23.6 and 572 ± 7.48 MPa, the UTSs are 676.7 ± 20.5 and 691 ± 15.9 MPa, and the elongations are 12.4 ± 2.1% and 17.9 ± 0.9% for C and S samples respectively, which has shown to be superior to the previous work with the same layer thickness [13]. Moreover, S samples performed better with smaller deviations in terms of strength and ductility in comparison with C samples.

To elucidate the reasoning for the different mechanical performance of C and S samples, μ-CT data and the corresponding cross-sectional OM images are shown in Fig. 2(a) to 2(b). With the same scanning parameters but a different scanning strategy, C sample contains more cracks/pores compared to S sample. Moreover, for S sample, the cracks are randomly distributed which is a common phenomenon for metal AM fabricated parts [32].

However, for C sample, the cracks align themselves along the build direction. This phenomenon is validated in the OM images with the aligned cracks are circled in red. It is worth noting that the chimney-like cracks were only present under the chessboard scanning strategy with stripe sizes of 1.25 and 2.5 mm. It does not appear under the stripe scanning strategy (stripe sizes from 0.5 mm to 10 mm) nor the chessboard scanning strategy with a stripe size of 5 mm. For tensile fracture surfaces shown in Fig. 2(c) and 2(d), more elongated cracks are detected in C sample and they are believed to cause its property worsening and deviation. Hot tearing/solidification cracks are observed in both samples (Fig. 2

(e) and 2(f))

showing typical protrusions with smooth surfaces [19].

It is reported that chemical heterogeneity at the high angle grain boundary helps to maintain the liquid film at the end of solidification, which causes solidification cracking of metal AM parts [19]. Elemental distribution in 3D form at atomic level was specifically conducted by means of APT technique near an intergranular crack as shown in Fig. 3(a). The APT tip specimen with a radius of ~50 nm was analysed (Fig. 3(b)). Fig. 3(c) reveals that distribution of the 4 constituting elements is uniform and their atomic ratios remain to be consistent. Obviously, no elemental segregation is seen at the intergranular crack. Thus, it is suggested that the solidification cracking could not be attributed to elemental segregation for this HEA. Fig. 4(a) schematically illustrates a hot tearing mechanism in SLM process. As the liquid metal solidifies after the laser melting for each layer, cellular structures grow towards different directions within each grain. Residual stress ascribed to thermal shrinkage acts as a tension force perpendicular to the grain boundary. The unsolidified semi-solid metal tends to fill up any open channel via gravity/capillary force, known as liquid feeding [33]. It is known that cracking susceptibility is related to the grain size for welding [34]. During solidification, grain boundaries are the weak "linkages" within the solid due to incoherent atomic packing.

Therefore, under a uniform tensile force, fine-grained microstructure provides more grain boundaries to share the load and cracking is less likely to occur. Given the large grain size of SLM-built CoCrFeNi, there are fewer grain boundaries to withstand the residual stress. To assess the residual stress in the as-built CoCrFeNi samples, TEM studies were carried out. To quantify the effect of grain size towards the hot tearing for the SLM processing of HEA, the Rappaz-Drezet-Gremaud (RDG) criterion is employed as it has been proven to be one of the most effective models in cracking susceptibility prediction [35]. The original RDG model considers the depression pressure within the mushy zone induced by both solidification shrinkage ( ) and mechanical strain ( ) [36]: with whereas is the thermal gradient, is the shrinkage factor, is the dynamic viscosity, is the characteristic microstructural length, taken as the cellular structure width here, is the solid-liquid front growth velocity, is the solid fraction, is the strain rate, is the elongation and is the liquidus temperature. is the coherency temperature, below this value, liquid feeding problem becomes trivial as liquid film breaks down into isolated drops [36]. Though the RDG model was built without taking the grain boundary effects into consideration, Wang et al. successfully used it to explain the effect of grain boundary angle towards the laser-welding cracking susceptibility [37]. Therefore, the present study aims to establish a new set of characteristic depression pressure values specifically for metal AM.

Unlike casting and welding, SLM involves an extremely high thermal gradient of ~10 6 K/m [26]. It makes the characteristic depression pressure significantly higher compared to the conventional values. Like the original calculations made by the RDG model, the strain rate is assumed to be a constant value of 10 -2 s -1 [36,[START_REF] Hagenlocher | Science and Technology of Welding and Joining[END_REF]. Due to the different nucleation tendencies and existing experimental data, it is assumed that for a cross-section normal to the build direction, the largest average grain size can be obtained for SLM-built CoCrFeNi, SS316L and Al7075 are 10 6 , 10 6 , and 10 4 μm 2 respectively [START_REF] Sistiaga | [END_REF]40]. To incorporate the effect of grain size, the strain rate is proportionally decreased as the lateral grain size decreases. Applying existing solidification path data (Fig. 4(c)) and their respective materials properties (Supplementary Table 1), we plot the depression pressure curves for CoCrFeNi, SS316L and Al7075, shown in Fig. 4(d). The coherency temperature is taken as the temperature when = 0.93 as suggested by [41]. Typically, is inversely proportional to solidification temperature range. However, even with a short solidification range of ~20°C for CoCrFeNi, it still has a relatively high depression pressure due to its high viscosity value of 3.15×10 -2 kgm -1 s -1 [START_REF] Mohr | TMS2019 Proccedings[END_REF]. It is about 1 order of magnitude higher compared to the common alloys (Supplementary Table 1). It is observed that as the grain size decreases, the contribution of mechanical strain decreases exponentially. By substituting in the existing experimental data, it is found that once the corresponding grain size exceeds a characteristic depression pressure value, solidification cracking will start to occur. It is worth noting that, when the grain size is close to the transition point, only minor cracking appeared (Supplementary Fig. 2(a)). Based on the above-mentioned results, it is conceived that hot tearing for the SLM-built CoCrFeNi could be overcome when its average grain size is below 10 4 μm 2 . Therefore, the focus of our future work will be refining the microstructure so as to conquer the hot tearing issue for AM processing of various HEAs.

In summary, this work revealed the origin of hot tearing for a HEA CoCrFeNi additively manufactured by SLM process. A wide range of processing parameters were implemented on pre-alloyed powder, but intergranular solidification cracks were present throughout the builds. The 4 constituting elements were found to be uniformly distributed on the periphery of cracks, which means that elemental segregation induced thin film of liquid may not be the reason for the hot tearing of this alloy. It is found out that severe residual stress due to the formation of large grains is the root cause and there will be a characteristic depression pressure limit beyond which solidification cracking tends to occur.

Fig. 4 (

 4 Fig. 4(b) shows a cellular structure in the <1-10> zone axis. The cell boundaries are marked by white arrows. Previously reported SLM-built cellular structures typically have a clean interior surrounded by dislocation cell boundaries [2]. However, the present as-built CoCrFeNi sample has dense entangled dislocations within the cell. It suggests that the SLMbuilt CoCrFeNi is under a high level of residual stresses.

Figure 1 .

 1 Figure 1. OM images of SLM-built HEA CoCrFeNi with (a) chessboard (denoted as C) and (b) stripe (denoted as S) scanning strategies. EBSD IPF color maps with respect to the build direction of the selected areas for (c) C sample and (d) S sample. Lack-of-fusion cracks and spherical pores are highlighted in black and red circles, respectively. Schematics of (e) chessboard and (f) stripe scanning strategies for SLM process.

Figure 2 .

 2 Figure 2. μ-CT images with their corresponding cross-sectional OM images of the (a) C and (b) S samples. The chimney-like crack in C sample is circled in red dashed line. The tensile fracture surfaces of (c) C and (d) S samples with solidification cracks highlighted by white arrows. Enlarged view of solidification cracking for (e) C and (f) S samples showing typical smooth protrutions.

Figure 3 .

 3 Figure 3. (a) SEM image showing that APT specimens were acquired on the periphery of a solidification crack. (b) SEM image showing an APT specimen tip with a dimension of ~50 nm. (c) APT reconstruction volumes showing atom distribution of the 4 constituting elements (Co, Cr, Fe and Ni).

Figure 4 .

 4 Figure 4. (a) A generic schematic illustration of the solidification cracking mechanism in AM process. (b) TEM image showing a large number of dislocations within the as-built cellular microstructure. Cell boundaries are indicated by white arrows. A selected area diffraction (SAED) pattern is shown with a <1-10> zone axis. (c) Solidification paths of SS316L, CoCrFeNi and Al7075. (d) Hot tearing propensity with respect to their grain size and typical depression pressures for the chosen SLM-built alloys.

Figure

  Figure(s) Click here to download high resolution image

Acknowledgements

The authors are grateful for the use of resources in Singapore Centre for 3D Printing funded by the National Research Foundation (NRF), Singapore.

Supplementary Material

Click here to download Supplementary Material: Supplementary material.docx


Hugo Cayuela

Jean-François Lemaître

Lorenzo Rugiero

Massimo Capula

Luca Luiselli

Asynchrony of actuarial and reproductive senescence: a lesson from an indeterminate grower

Aging evolutionary theories predict that patterns of actuarial and reproductive senescence should be aligned, with a common onset of senescence set at the age of first reproduction. However, a few empirical studies reported asynchrony between actuarial and reproductive senescence. This asynchrony is expected to be particularly pronounced in organisms with indeterminate growth. Yet, this process is still poorly documented due to the lack of long-term demographic data on knownaged individuals. We investigated the asynchrony of actuarial and reproductive senescence in the

Introduction

Although a few studies have challenged the idea of an inevitable senescence in nature (e.g. Martinez 1998, Cayuela et al. 2020), it is now broadly recognized that senescence, the decline in the contribution to fitness with increasing age, is a widespread phenomenon across the animal kingdom (Nussey et al. 2013;[START_REF] Gaillard | An integrative view of senescence in nature[END_REF]. Senescence can be due to an age-specific decrease in survival probabilities (i.e., actuarial senescence) and/or reproductive performance (i.e., reproductive senescence), two processes usually anchored in the same theoretical framework [START_REF] Gaillard | The Williams' legacy: A critical reappraisal of his nine predictions about the evolution of senescence[END_REF]. Evolutionary theories predict that patterns of actuarial and reproductive senescence should be aligned (Williams 1957, Maynard-Smith 1962), with a common onset of senescence at the age of first reproduction (Williams 1957). However, few recent studies, performed in both laboratory and free-ranging conditions, reported asynchrony between actuarial and reproductive senescence patterns (reviewed in [START_REF] Gaillard | The Williams' legacy: A critical reappraisal of his nine predictions about the evolution of senescence[END_REF]. The factors modulating the degree of synchrony are still poorly understood.

The asynchrony between actuarial and reproductive senescence could be particularly pronounced in organisms with indeterminate (i.e., continuous) growth (e.g., fishes, amphibians, and reptiles) where survival and reproductive traits (e.g., fecundity, offspring survival) are often positively correlated with maternal body size, and thus indirectly determined by female age [START_REF] Hoekstra | The untapped potential of reptile biodiversity for understanding how and why animals age[END_REF]. This asynchrony could result from non-linear relationships between fitness components and age. The fecundity of females with indeterminate growth recurrently increases with body size and age in a linear way (e.g., in reptiles, Jenkins et al. 2009[START_REF] Setser | Reproductive ecology of female Mexican lance-headed rattlesnakes[END_REF][START_REF] Rose | Reproductive frequency and size-dependence of fecundity in the giant gartersnake (Thamnophis gigas)[END_REF], whereas the strength and the form of the relationship between mortality and age after sexual maturity can be highly variable among species (Colchero et al. 2019, Cayuela et al. 2019a) and populations (Cayuela et al. 2020).

To date, the asynchrony of actuarial and reproductive senescence in indeterminate growers is still poorly documented, mainly due to the lack of long-term demographic data on known-aged day lasting from 0800 to 1800), by two or three people who searched intensively for snakes in the dilapidated building and neighbouring area of about 100× 80 m. Snakes were captured by hand, and individually marked by ventral scale clipping for future identification. At their first capture, female snakes were measured for snout-vent-length (SVL). We did not measure them on further captures to minimize disturbance and hence risks of their departure from the nesting area, as permanent emigration would bias capture-recapture inferences. We were able to assess the age of 16 females that were captured at juvenile stage over the capture-recapture survey. Detailed information about the number of captures performed each year is given in Table S1.

Several females (N = 16) on their first capture were palpated in the abdomen to verify their pregnancy status and to count the number of eggs [START_REF] Filippi | Long-term fidelity to communal oviposition sites in Hierophis viridiflavus[END_REF]). We also captured 44 additional females in the surroundings of the study area to increase the number of observations of egg clutch size (N = 60 in total).

Reproductive senescence: assessing age-dependent variation in clutch size

First, we verified that SVL was a good proxy of female age by quantifying the correlation between age and SVL in a small set of females for which the age was known (N = 16). Then, we examined how annual egg clutch size was associated with SVL in a larger sample of females (N = 60). In both cases, we used regression models where SVL or egg clutch size were included as the response variable and body size or age were included as explanatory variables. We compared null (i.e. without covariate), linear, and polynomial models using their AIC. We graphically verified the fit of the model by examining model residuals. The analyses were conducted in R program (R Development Core Team 2018).

Modeling age-dependent mortality using BaSTA

Reproductive senescence pattern

Our results confirmed that female SVL linearly increased with age (Fig. 1A and Table S2 for model selection procedure) and was thus a robust proxy of their age (! " = 0.43). Next, we showed that annual egg clutch size increased with SVL (Fig. 1B). Although the polynomial model was close (less than 2 AIC points between the two models), the linear model was selected ! " (0.74) (Table S3). Overall, our analyses highlighted that female annual fecundity increased linearly with SVL and age.

Actuarial senescence pattern

BaSTA models revealed that capture-recapture data were best described by a logistic function without any refinement of the shape (Table 1). The cumulative probability of surviving was 0.75 until age six, 0.50 until age seven, and 0.25 until age eight (Fig. 1C). Furthermore, the model showed a strong positive effect of age on hazard rate (Fig. 1D), which indicates an early and strong actuarial senescence in H. viridiflavus.

Discussion

Using capture-recapture data collected over a 29-year period, we showed that H. viridiflavus female did not experience any fecundity loss late in life, suggesting a negligible reproductive senescence.

In contrast, they suffered from an early and severe actuarial senescence. Taken together, those results revealed a pronounced asynchrony in the senescence process of two main fitness components.

general inferences about the ecological and biological factors modulating the asynchrony between aging components. At the intraspecific level, abiotic factors such as temperature and habitat predictability regulate lifespan and actuarial senescence [START_REF] Sears | Geographic variation in the life history of the sagebrush lizard: the role of thermal constraints on activity[END_REF][START_REF] Hjernquist | Seasonality determines patterns of growth and age structure over a geographic gradient in an ectothermic vertebrate[END_REF], Zhang et al. 2018, Cayuela et al. 2019b), as well as reproductive investment (Morrison et al. 2003, Cayuela et al. 2016[START_REF] Hughes | Latitudinal variation in life history reveals a reproductive advantage in the Texas Horned Lizard (Phrynosoma cornutum)[END_REF], in ectotherm vertebrates with an indeterminate growth, which could ultimately modify the asynchrony of senescence patterns. In our study system, the low abundance of natural predators compared to other habitats occupied by H. viridiflavus (Philippi et al. 2007) could also buffer actuarial senescence by reducing predation-induced mortality late in life, which might lead to intraspecific variation in the degree of asynchrony between aging components. Furthermore, covariation between mortality patterns and several phenotypic traits of indeterminate growers could result in interspecific variation in the asynchrony of actuarial and reproductive senescence. In particular, the magnitude of the asynchrony could be diminished in species with size-dependent fecundity and negligible or even "negative" actuarial senescence (sensu [START_REF] Vaupel | The case for negative senescence[END_REF], Jones & Vaupel 2017) Relationship between hazard rate and age described by a logistic function.

  allowed by high regeneration capacities (e.g., salamanders, Cayuela et al. 2019a), or chemical (e.g., venomous snakes, Cayuela et al. 2019a, Colchero et al. 2019) and armour protection (e.g., turtles, Warner et al. 2016). Overall, our study raises important questions about the ecological and evolutionary mechanisms allowing the uncoupling of reproductive and actuarial senescence processes.

Fig. 1 .

 1 Fig.1. Reproductive and actuarial senescence in females of Hierophis viridiflavus. (A) Linear

Acknowledgments

We thank Christopher R. Friesen and two other anonymous reviewers for their helpful comments.


Hugo Cayuela
email: hugo.cayuela51@gmail.com

Jean-François Lemaître

Eric Bonnaire

Julian Pichenot

Benedikt R Schmidt

Population position along the fast-slow life-history continuum predicts intraspecific variation in actuarial senescence

Introduction

Most commonly accepted evolutionary theories of ageing posit that survival should decline with increasing age in any age-structured population (Hamilton 1966), a demographic process known as actuarial senescence. [START_REF] Medawar | An unsolved problem of biology[END_REF] mutation accumulation theory first stated that organisms age because the strength of natural selection weakens with age after first reproduction and therefore there is no purging of deleterious mutations that are only expressed late in life. In addition, actuarial senescence can emerge as a by-product of natural selection through antagonistic pleiotropy [START_REF] Williams | Pleiotropy, natural selection, and the evolution of senescence[END_REF]). An allele may confer a benefit to the bearer early in life but may also be responsible for an impaired survival later in life. Finally, the disposable soma theory of aging postulates that actuarial senescence can result from a trade-off between allocation to reproduction in early life and somatic maintenance (Kirkwood 1977, Kirkwood & Austad 2000). In other words, individuals that preferentially allocate resources to growth and/or reproduction (e.g., gamete production and parental care) will have much less resources for somatic maintenance (e.g., enzymebased repair mechanisms), which will ultimately lead to a decline in performance of fitness-related traits (e.g., survival) at advanced ages. So far, most studies focused on aging in the wild have been embedded within the two theoretical frameworks offered by both antagonistic pleiotropy and disposable soma theories of aging as these two theories share the similar prediction of a trade-off between reproductive effort and actuarial senescence [START_REF] Gaillard | The Williams' legacy: A critical reappraisal of his nine predictions about the evolution of senescence[END_REF].

Although a few case studies on free-ranging populations failed to detect any increase in mortality rate with age (Jones et al. 2014, Cayuela et al. 2019c), most recent syntheses revealed two major facts. First, actuarial senescence is a nearly ubiquitous process in the living world [START_REF] Nussey | Senescence in natural populations of animals: widespread evidence and its implications for bio-gerontology[END_REF], Shefferson et al. 2017). Second, patterns of senescence can be highly variable among species (Jones et al. 2014[START_REF] Tidière | Comparative analyses of longevity and senescence reveal variable survival benefits of living in zoos across mammals[END_REF][START_REF] Colchero | The diversity of population responses to environmental change[END_REF]. Comparative analyses have shown that senescence patterns across multi-cellular organisms can be predicted by ecological traits, lifestyles and covariation among life-history traits [START_REF] Péron | Age at the onset of senescence in birds and mammals is predicted by early-life performance[END_REF][START_REF] Ricklefs | Life-history connections to rates of aging in terrestrial vertebrates[END_REF], Gaillard et al. 2016, Salguero-Gómez & Jones 2017). In particular, both age at the onset of senescence and rates of senescence appear to be linked to the position of a species along the fast-slow life-history continuum. Organisms that occupy the fast end of the continuum -short generation time, high annual fecundity, and low mean survival rates (Stearns 1983) -tend to experience earlier and faster actuarial senescence than organisms at the slow end of the continuum (Jones et al. 2008a, Salguero-Gómez & Jones 2017). Since life history variation does not only occur among, but also within species [START_REF] Berven | Interpreting geographic variation in life history traits[END_REF], Cayuela et al. 2017), intraspecific variation in senescence patterns is expected and can be selected for (Stearns 2000, Stearns et al. 2000, Brommer et al. 2007, Holand senescence and the speed and the shape of the relationship between mortality rate and age) differed between populations with slow and fast life histories. We focused our analyses on postmetamorphic survival as we did not expect senescence during larval development and assumed that senescence can begin once structural development is completed (i.e., metamorphosis). We expected that the pattern of intraspecific variation in actuarial senescence among populations would be similar to pattern in variation among species. Therefore, we tested the prediction that a fast paceof-life should be associated with an earlier, faster actuarial senescence among populations of B. variegata.

Material and methods

Study system

Bombina variegata populations occur in different types of habitats where the spatiotemporal pattern of breeding resource availability differs widely. In riverine environments, the patches of rock pools used by toads to reproduce are constantly available in space and time, making breeding resources highly predictable at the scale of the lifetime of an individual. By contrast, rut patches resulting from logging operations appear and disappear stochastically in forests, making breeding resource availability unpredictable. The annual probability of patch appearance varies from 0.20 to 0.50 while the rate disappearance ranges from 0.05 to 0.20, depending on the year and the population (Cayuela et al. 2016a).

Previous studies showed that the populations from these two types of environment have different life history strategies. Fast populations from forest environment have accelerated life histories, with lower postmetamorphic survival and higher fecundity than in riverine habitats.

Moreover, the fast life history is associated with higher dispersal (i.e., dispersal syndrome) caused by the high patch turn-over resulting from logging operations (Cayuela et al. 2016c(Cayuela et al. , 2019b)).

Dispersal is facilitated by a phenotypic specialization that is transgenerationally transmitted: individuals from fast forest populations reared in common garden have a higher boldness and exploration propensity than those from riverine slow populations (Cayuela et al. 2019b). The association between high dispersal maintained throughout lifespan, pre-natally determined behavior, and accelerated life history (i.e., life history syndrome; [START_REF] Réale | Personality and the emergence of the pace-of-life syndrome concept at the population level[END_REF]) is likely determined by polygenic bases and gene with pleiotropic effects (Saastamoinen et al. 2018). Furthermore, changes in dispersal regime within spatially structured populations of B. variegata lead to strong variation in neutral genetic structure. Fast forest populations have a higher genetic pattern used to identify individuals is not fixed before that age. Yet, this method permits to survey the toads for the rest of their lives in a non-invasive way. Only individuals can be aged which are one or two years old at first capture (Cayuela et al. 2016b), leading to partially known birth date; this incomplete knowledge of age of individuals is handled using Bayesian survival trajectory analyses (see the section Capture-recapture modeling). Sex cannot be assessed with certainty before sexual maturity (2-3 years old) due to the lack of nuptial pads in immature males. For this reason, we did not include an effect of sex in the analyses. Yet, we expect that excluding sex from our analyses should not alter our conclusions since previous studies in these populations did not detect sex-specific effects on survival (Cayuela et al. 2016a(Cayuela et al. , 2016c)). Details about the number of survey years, individual age and sampling effort are provided in Appendix 1 (Table 2). A more complete description about the capture-recapture survey, the individual recognition method and population description can be found in Cayuela et al. (2016aCayuela et al. ( , 2016c)).

Capture-recapture modeling

We investigated actuarial senescence patterns using Bayesian survival trajectory analyses implemented in the R package BaSTA (Colchero et al. 2012a(Colchero et al. , 2012b)). BaSTA allowed us to account for imperfect detection, left-truncated (i.e., unknown birth date (age)) and right-censored (i.e., unknown death date) capture-recapture data in our analysis. Our analyses focus on the postmetamorphic stage at which senescence is expected to occur (as in [START_REF] Colchero | The diversity of population responses to environmental change[END_REF]. BaSTA allows estimation of two parameters: age-dependent survival and the proportion of individuals dying at a given age (i.e. age-dependent mortality rate). Given the results of previous analyses (Cayuela et al. 2016a(Cayuela et al. , 2016b)), we allowed recapture probabilities to vary among years. As the study period and number of survey years differ among populations (Appendix 1), the four populations were analyzed separately. We used deviance information criterion (DIC) to select models that fitted the data best [START_REF] Colchero | The diversity of population responses to environmental change[END_REF]) and we compared the outputs of the best-supported model of the four populations by inspecting mean estimates and 95% credible intervals (CRI) [START_REF] Anderson | Suggestions for presenting the results of data analyses[END_REF][START_REF] Amrhein | Retire statistical significance[END_REF]. This allowed us to investigate population-specific variation in the shape of the age-specific mortality patterns. We considered the four mortality functions implemented in BaSTA: exponential, Gompertz, Weibull and logistic. For the three last functions, we considered three potential shapes: simple that only uses the basic functions described above;

Makeham [START_REF] Pletcher | Model fitting and hypothesis testing for age-specific mortality data[END_REF]; and bathtub (Silver 1979). As individuals which are less than one year old cannot be individually recognized, we conditioned the analyses at a minimum age of one. Four MCMC chains were run with 50000 iterations and a burn-in of 5000. Chains were thinned by a after individuals are 2-3 years old and reached an asymptote at an age of 3 to 4 years. This suggests an early, fast senescence with a short reproductive lifespan.

Age-dependent survival patterns were relatively similar between the populations with the same life history strategies but differed markedly between strategies (Fig. 3). In slow populations (Fig. 3A and3B), the cumulative probability of surviving until a given age decreased slowly over a toad's lifetime: it was 0.78 until age three, 0.43 until age six, 0.17 until age nine, and finally 0.00 until age 12 (Appendix 1, Table 23). In fast populations (Fig. 3C and3D), cumulative survival probability decreased rapidly after two years: it was 0.58 until age three, 0.07 until age six, and 0.00 until age nine (Appendix 1, Table 45.

Discussion

We provide the first clear evidence for actuarial senescence in a wild amphibian that is associated with an intraspecific slow-fast life history difference. The onset of actuarial senescence was earlier in individuals from fast populations and individuals show a faster rate of actuarial senescence than individuals in slow populations.

Population viability, elasticity of demographic parameters, and potential confounding factors

Population matrix models showed that the four populations had population growth rates higher than 1 and can therefore be considered to be viable. Moreover, we highlighted changes in the elasticity of demographic rates: slow population have higher elasticities for adult survival and lower elasticities for pre-breeding survival and fecundity than fast populations. Taken together, our results suggest that life history acceleration in forest populations is not caused by an ongoing demographic decline or transient dynamics. By contrast, variation in life history speed seems to rely on stable demographic strategies covarying with behavioral syndromes that likely have genetic bases (Cayuela et al. 2019b).

Although the habitat type (forest vs riverine) is confounded with latitude in our study system, it is unlikely that the difference in latitude has an effect on demographic traits and actuarial senescence. First, a previous study showed that climate has idiosyncratic effects on population dynamics and there was no evidence for an influence of latitude on demographic rates in B. variegata (Cayuela et al. 2016b). Second, populations with slow life history strategies occurring in quarry habitats can be found in the northern part of the study area (Cayuela et al. 2019a), which suggests that latitude is not an important predictor of demographic rates. Third, life history variation precise quantification of this onset) and had a faster rate of actuarial senescence than individuals from slow populations. In our study system, individuals from fast populations occurring in forest environment face higher level of environmentally-driven mortality (associated with the unpredictability of breeding patches) than individuals from the slow populations. In fast B. variegata populations, logging management frequently lead to the destruction of breeding patches, which decreases adult survival via direct mortality during harvesting operations (Cayuela et al. 2018). In addition, as the natural silting of ruts decreases local breeding success over time, individuals frequently disperse to colonize new patches and likely incur energy and mortality costs related to the different stages of the dispersal process (i.e., emigration, transience, and immigration;Boualit et al. 2019, Cayuela et al. 2019b). By contrast, slow populations occur in non-managed habitats with very limited human activity where dispersal is virtually absent, which likely decrease environmentally-driven mortality (Cayuela et al. 2019b). Since [START_REF] Williams | Pleiotropy, natural selection, and the evolution of senescence[END_REF] pioneering contribution on this topic, the relationship between environmentally-driven mortality and actuarial senescence has been heavily discussed (see Gaillard &[START_REF] Gaillard | The Williams' legacy: A critical reappraisal of his nine predictions about the evolution of senescence[END_REF]Ronget et al. 2017 for reviews). According to [START_REF] Williams | Pleiotropy, natural selection, and the evolution of senescence[END_REF], high level of adult mortality should lead to a much stronger rate of actuarial senescence, a prediction which was validated by simulations [START_REF] Gaillard | The Williams' legacy: A critical reappraisal of his nine predictions about the evolution of senescence[END_REF], at least in absence of densitydependent or condition-dependent mortality over the entire life course [START_REF] Abrams | Does increased mortality favor the evolution of more rapid senescence[END_REF][START_REF] Williams | Antagonistic pleiotropy, mortality source interactions, and the evolutionary theory of senescence[END_REF]; see also [START_REF] Moorad | Evolutionary Ecology of Senescence and a Reassessment of Williams' 'Extrinsic Mortality' Hypothesis[END_REF] for a recent review). For instance, if mortality during the juvenile period is strongly condition-dependent, only the most robust individuals will reach adulthood which can postpone and/or decrease actuarial senescence rates (Ronget et al. 2017). Such information is currently lacking in B. variegata but would provide important insights on how environmentally-driven mortality modulate the selective forces acting against actuarial senescence. In our study system, actuarial senescence likely evolves from an adjustment of resource allocation strategies in response to environmentally-driven mortality (Kirkwood & Rose 1991;Stearns 2000[START_REF] Baudisch | Getting to the root of aging[END_REF]. In others words, individuals from the fast populations might be able to cope with higher level of environmentally-driven mortality by allocating more resources to reproduction at the expense of somatic maintenance (see Cayuela et al. 2016aCayuela et al. , 2016b for an in-depth discussion). This might ultimately lead to a faster decline in survival probabilities of individuals with age, as expected under the disposable soma theory of ageing (Kirkwood & Austad 2000). In contrast, environmentally-driven mortality is lower in riverine environment where the predictability of breeding patches is high and human activity is virtually inexistent.

Accordingly, individuals may allocate fewer resources to reproduction -females produce a lower number of juveniles per year (Cayuela et al. 2016a) -and more to somatic maintenance, resulting Table 2. Model selection results. Deviance information criterion (DIC) for each of the mortality function considered in the four studied populations of Bombina variegata. We considered the four mortality functions implemented in BaSTA program: exponential (EXP), Gompertz (GOM), Weibull (WEI) and logistic (LOG). For the three last functions, we considered three potential shapes: simple that only uses the basic functions described above ("simple"); Makeham ("make"); and bathtub ("bath").

Function

Fig. 2 .

 2 Fig. 2. Mortality rate (i.e., proportion of individuals dying at a given age) in two fast (FAST1 and FAST2) and slow (SLOW1 and SLOW2) populations of Bombina variegata. The predictions on the left of the vertical dashed line correspond to observed ages while those on the right are model projections.

Fig. 3 .

 3 Fig. 3. Survival probability until a given age in two fast (FAST1 and FAST2) and slow (SLOW1 and SLOW2) populations of Bombina variegata. The predictions on the left of the vertical dashed line correspond to observed ages while those on the right are model projections.

Acknowledgments

We thank all the field workers that have contributed to data collection and Erin Muths for her comments on the manuscript. This research was funded by the Lorraine DREAL, the Rhône-Alpes DREAL, the Agence de l'Eau Rhône-Alpes, the Agence de l'Eau Rhin-Meuse, the Office National des Forêts, the Conseil Régional de Lorraine, the Conseil Régional de Champagne-Ardenne, the Conseil Régional de Picardie, the Conseil Général de l'Aisne, the Conseil Général d'Ardèche, the Conseil Général d'Isère and the Communauté de Communes de l'Argonne Ardennaise (2C2A).

Toad capture was authorized by the Préfecture de l'Ardèche (arrêté no. 2014-288-002) and the Préfecture de la Meuse (arrêté no. 2008-2150).

2017), then phylogenetic constraints would be a weak explanation for senescence patterns across the tree of life [START_REF] Antonovics | Ontoecogenophyloconstraints? The chaos of constraint terminology[END_REF]. Second, if there is intraspecific variation in senescence, this variation could explain the weakness of the phylogenetic signals in actuarial senescence patterns in several taxa (Salguero-Gómez & Jones 2017). It may well be that phylogenetic signals could become more apparent once intraspecific variation is taken into account.

Such a combination of microevolutionary and macroevolutionary patterns of senescence would lead to a deeper understanding of the evolutionary biology of senescence.


Zheng Zhang
email: zhangzheng@aircas.ac.cn

Ping Tang
email: tangping@aircas.ac.cn

Thomas Corpetti
email: thomas.corpetti@univ-rennes2.fr

Time Adaptive Optimal Transport : A Framework of Time Series Similarity Measure

Keywords: Optimal transport, time series, similarity measure, Sinkhorn distance, classification

Similarity measure is a critical tool for time series analysis. However, currently established methods, for instance, dynamic time warping (DTW) and its variants, are still facing some issues such as non-maximum-to-maximum alignment and pathological alignment, etc. Despite many attempts to improve, these issues remain stubborn because they are directly caused by the intrinsic mechanism of DTW. Thinking out of the context of DTW based methods, we propose in this paper a new time series similarity measure framework which we call Time Adaptive Optimal Transport (TAOT). As its name implies, TAOT is based on optimal transport, a powerful distance measure for histograms and probability distributions, and TAOT inherits several promising properties from optimal transport to tackle the problems of classic DTW based methods. We make optimal transport capable of handling time series data by considering both observed values and their corresponding time coordinates simultaneously. TAOT can generate a many-to-many alignment between time series that further releases the search space for a more correct result. Experimental results show that TAOT can outperform other widely used similarity measures on classification tasks on multiple datasets. We also introduce the parameter extracting and visualization strategies of TAOT in this paper.

I. INTRODUCTION

S IMILARITY measure is fundamental to many machine learning and data mining tasks, such as classification [START_REF] Chen | Similarity-based classification: Concepts and algorithms[END_REF]- [START_REF] Wan | Study on the satellite telemetry data classification based on self-learning[END_REF], clustering [START_REF] Liao | Clustering of time series data ąła survey[END_REF]- [START_REF] Li | Time-series representation and clustering approaches for sharing bike usage mining[END_REF], indexing [START_REF] Bartolini | Warp: Accurate retrieval of shapes using phase of fourier descriptors and time warping distance[END_REF]- [START_REF] Li | Speed up similarity search of time series under dynamic time warping[END_REF], etc. In the context of time series analytics, similarity measure has long been a research hotspot and many methods have been proposed [START_REF] Lhermitte | A comparison of time series similarity measures for classification and change detection of ecosystem dynamics[END_REF], [START_REF] Wang | Experimental comparison of representation methods and distance measures for time series data[END_REF]- [START_REF] Savaş | Evaluation of different algorithms for measuring the similarities of trajectory datasets[END_REF]. These similarity measures can be roughly categorized into: time-rigid measures (Euclidean distance), time-flexible measures (dynamic time warping) [START_REF] Folgado | Time alignment measurement for time series[END_REF], [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF], feature-based measures (Fourier coefficients) [START_REF] Oppenheim | Discrete-time signal processing[END_REF]- [START_REF] Zhang | Local feature based dynamic time warping[END_REF], and model-based measures (auto-regression and moving average model) [START_REF] Serri | An empirical evaluation of similarity measures for time series classification[END_REF], [START_REF] Fu | A review on time series data mining[END_REF]. Among these methods, dynamic time warping (DTW) [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] and its variants [START_REF] Li | Similarity measure based on incremental warping window for time series data mining[END_REF]- [START_REF] Zhang | Dynamic time warping under limited warping path length[END_REF] are probably the most popular and established ones [START_REF] Wang | Experimental comparison of representation methods and distance measures for time series data[END_REF], [START_REF] Fu | A review on time series data mining[END_REF], [START_REF] Ding | Querying and mining of time series data: experimental comparison of representations and distance measures[END_REF]. However, this kind of DTW-based methods still exhibit several intrinsic alignment problems: a) Pathological alignment: A critical job for similarity measure of time series is to tackle time distortions. DTW screens most time distortions by realigning data points in the two time series to be compared, and that helps DTW outperform Euclidean distance in many scenarios [START_REF] Ding | Querying and mining of time series data: experimental comparison of representations and distance measures[END_REF]. However, the intrinsic rules used by DTW to generate the new alignment can lead to pathological alignment, where a single point in one time series links to a large subsection of another time series [START_REF] Keogh | Derivative dynamic time warping[END_REF], as shown in Fig. 1(a). In order to avoid this undesirable alignment, various constraints for DTW [START_REF] Itakura | Minimum prediction residual principle applied to speech recognition[END_REF]- [START_REF] Berndt | Using dynamic time warping to find patterns in time series[END_REF] have been proposed, for instance, Sakoe-Chiba band and Itakura parallelogram. But most of these constraints are rigid such that they take a risk of preventing the correct alignment from being generated. b) Pairing of maximum values: In some applications, the maximum value of a time series is of overwhelming importance, and thus it is decisive to guarantee that the maximum value of one time series aligns to the maximum value of another time series. Unfortunately, DTW and its existing variants provide no such guarantee, and to some extent this requirement even conflicts with the underlying mechanism of DTW. A mandatory pairing of maximum values can break the time-monotonicity of DTW and prevent DTW from finding the optimal alignment with the lowest global cost. Other common similarity measures, such as Euclidean distance and Minkowski distance, can only succeed in a certain situation when two peaks happen at the same time point, because they obey a time-rigid alignment rule. c) Many-to-many alignment: From the perspective of typology, there should be three categories of alignments between time series, one-to-one, one-to-many, and many-tomany, as illustrated in Fig. 1(b-c-d) respectively. Among the three categories, the potential of one-to-one alignment and one-to-many alignment has already been exploited; for instance, time-rigid methods employ one-to-one alignment and DTW-based methods employ one-to-many alignment. However, the use of many-to-many alignment is still an open direction. In theory, many-to-many alignment has a larger capacity and a higher flexibility. Therefore, many-to-many alignment based methods deserve to be developed. d) Privilege of time ordering: The similarity between time series is measured mainly from two aspects: observed values and time ordering. To some extent existing methods always consider time ordering in the first place. For instance, DTW first sets a rule of time ordering (continuity, monotonicity, etc), and then finds the minimum accumulated cost of observed values under the given time ordering rule. Euclidean distance simply follows a time-rigid ordering. This privilege of time ordering narrows the search space and as a result forfeits the opportunity to find a more correct alignment. Instead, a new method that considers observed values and time ordering simultaneously may further release the capacity and lead to a more accurate similarity measure.

A review of literature shows a lack of method effectively responding to the aforementioned issues, and one reason is that these alignment problems are closely related to the core mechanism of DTW. Thinking outside the context of classic DTW-based methods, we find Optimal Transport (OT) [START_REF] Villani | Optimal transport: old and new[END_REF], also known as the Earth Mover's Distance (EMD) [34] or Wasserstein Distance, is a successful method to compare probabilities or histograms [START_REF] Ling | An efficient earth mover's distance algorithm for robust histogram comparison[END_REF]- [START_REF] Courty | Optimal transport for domain adaptation[END_REF], and theoretically it has several promising properties [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] to solve the aforementioned alignment problems of DTW. For example, OT does not rely on the time order of data points and it always aligns maximum values with each other. In the past the utility of OT was greatly limited by its high computational complexity, but fortunately this issue has been largely alleviated by one of its variants, Sinkhorn distance [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], which transfers OT into a problem with a fast solution by adding an entropic regularization term and ends up making the computation several orders of magnitude faster. Another barrier to the use of OT on time series data is that OT only considers the observed values while ignores the corresponding time coordinates.

In this paper, we propose a time-adaptive version of OT (TAOT) to serve as a new similarity measure framework between time series. We make Sinkhorn distance capable of handling time series data by considering both observed values and their corresponding time coordinates simultaneously. Classification experiments on multiple real-world and synthetic datasets are conducted to show the accuracy of TAOT and how the aforementioned alignment problems of DTW are coped with. Other closely related topics such as the parameter extracting strategies and the visualization of TAOT will also be discussed. The rest of the paper is structured as follows. Section II briefly revisits related works on dynamic time warping (DTW), optimal transport (OT), and their major variants. Section III details the proposed new similarity measure, namely TAOT. Section IV evaluates the performance of TAOT with classification experiments on UCR time series datasets and discusses several practical issues. Finally, Section V summarizes our main conclusions.

II. BACKGROUND A. DYNAMIC TIME WARPING

DTW was first introduced in the field of speech recognition in the influential work by [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] and soon become widely used because it can cope with time distortions effectively. DTW aims at finding the optimal alignment between time series that can achieve the minimum accumulated cost. Given two time series, A = a 1 , a 2 , ..., a i , ..., a m and B = b 1 , b 2 , ..., b j , ..., b n , DTW first constructs an m-by-n cost matrix D. Each matrix cell (i, j) contains the distance d(i, j) between the two data points a i and b j . Squared Euclidean distance is normally used when calculating the cost matrix, and thus d(a i , b j ) = (a i -b j )

2 . Then DTW tries to find the optimal alignment that leads to the minimum accumulated cost. The alignment is represented by a warping path W = w 1 , w 2 , ..., w k , ..., w K that consists of a continuous set of matrix cells. Each matrix cell, for instance, w k = (i, j) corresponds to a pairing of points a i and b j , and the cost between them d(i, j) is added to the final accumulated cost. Fig. 2 shows an example of the cost matrix and warping path. Eq. 1 defines the abovedescribed DTW distance:

DTW satisfies three time ordering rules: boundary condition, continuity, and monotonicity [START_REF] Keogh | Derivative dynamic time warping[END_REF]. The boundary condition restricts the warping path to start at the lower-left corner of the warping matrix w 1 = (1, 1) and finish at the upper-right corner w K = (m, n). The continuity restricts the allowable steps to adjacent cells. The monotonicity prevents all steps in the warping path from turning backward in any circumstances.

In practice, DTW distance can be efficiently calculated by the following recursive formula: dtw(i, j) = d(i, j) + min{dtw(i -1, j -1), dtw(i, j -1), dtw(i -1, j)}, where dtw(i, j) is the accumulated cost found in current cell (i, j) and the final DTW distance is dtw(m, n).

B. VARIANTS OF DTW

Despite its popularity, DTW was proposed decades ago and its pathological alignment problem has long been noted. In order to alleviate the problem, different variants of DTW have been proposed and we briefly review them here. These methods can be classified into two major categories.

The first category sets constraints on DTW. The most common and straightforward constraint is windowing [START_REF] Itakura | Minimum prediction residual principle applied to speech recognition[END_REF]- [START_REF] Berndt | Using dynamic time warping to find patterns in time series[END_REF], where allowable elements of warping path must be located within a warping window, as illustrated by the dashed line in Fig. 2. Different warping windows have been proposed, among which Sakoe-Chiba band and Itakura parallelogram are two most commonly used ones. In recent years, some researchers suggest to learn adaptive-shaped windows [START_REF] Yu | Dynamic time warping constraint learning for large margin nearest neighbor classification[END_REF], [START_REF] Candan | sdtw: computing dtw distances using locally relevant constraints based on salient feature alignments[END_REF] from the data instead of fixed-shaped windows. Another frequently used type of constraints is weighting. For example, Weighted DTW (WDTW) [START_REF] Jeong | Weighted dynamic time warping for time series classification[END_REF] is a representative method of weighting which applies different weights to temporally adjacent points when computing DTW. By penalizing further points, WDTW prevents minimum distance distortions caused by outliers and enhances the detection of similarity between two time series. Besides windowing and weighting, recently proposed LDTW [START_REF] Zhang | Dynamic time warping under limited warping path length[END_REF] constrains the maximum allowable length of the warping path based on the observation that pathological alignment always happens concurrently with an unusually long warping path. SP-DTW and its kernelization version SP-KrDTW [START_REF] Soheily-Khah | Sparsification of the alignment path search space in dynamic time warping[END_REF] address the sparsification of the alignment path search space for DTWlike measures to improve their efficiency without loosing accuracy.

The second category replaces the feature DTW considers. For example, Piecewise DTW (PDTW) [START_REF] Keogh | Scaling up dynamic time warping for datamining applications[END_REF], [START_REF] Cai | Piecewise statistic approximation based similarity measure for time series[END_REF] proposes to use a compact abstraction instead of the raw data to compute DTW, with the aim of avoiding the impact of outliers and accelerating the computation. PDTW first splits time series into fixed-length segments and compute the mean of each segment. Then these mean values are used instead of raw data points. A challenge of PDTW is the choice of the optimal size of segments. To avoid brute-force search, Parameter Free Piecewise DTW (FDTW) [START_REF] Siyou Fotso | Parameter Free Piecewise Dynamic Time Warping for time series classification[END_REF] proposes a heuristic search of the size of segments for PDTW on the basis of classification accuracy. Derivative DTW (DDTW) [START_REF] Keogh | Derivative dynamic time warping[END_REF] considers the local derivative of each data point rather than the raw value. Many variants of DTW can have a derivative version, for instance, WDTW and Weighted Derivative DTW (WDDTW) [START_REF] Jeong | Weighted dynamic time warping for time series classification[END_REF]. SC-DTW [START_REF] Zhang | Dynamic time warping under pointwise shape context[END_REF] views time series as 2D contours and employs shape context, a rich shape descriptor, to compute DTW. Local feature based DTW (LFDTW) [START_REF] Zhang | Local feature based dynamic time warping[END_REF] proposes a general framework to use different type of local features to generate the warping path.

C. OPTIMAL TRANSPORT

In the context of machine learning, OT [START_REF] Villani | Optimal transport: old and new[END_REF], [34](also known as Earth Mover's Distance (EMD) or Wasserstein Distance) has long been a powerful tool to compare probabilities or histograms [START_REF] Courty | Optimal transport for domain adaptation[END_REF], [START_REF] Kolouri | Optimal mass transport: Signal processing and machine-learning applications[END_REF], [START_REF] Robin | Detecting changes in forced climate attractors with wasserstein distance[END_REF]. OT is modeled as the solution to the transportation problem. Suppose there are a collection of mines mining iron ores, and a collection of factories that consume the iron ores. Given the amount of supply and demand of each mine and factory, and the shipment cost from each mine to each factory, OT can find the optimal allocation plan with a minimum total shipment cost for resolving the VOLUME 4, 2016 FIGURE 3. An example of optimal transport. We can observe that OT generates a many-to-many alignment.

supply-demand transports.

Given two probability distributions denoted as:

(A|p A ) = {(a 1 |p a1 ), (a 2 |p a2 ), ..., (a d |p ad )} (B|p B ) = {(b 1 |p b1 ), (b 2 |p b2 ), ..., (b d |p bd )} (2)
where a i or b i is the i-th observed value in the respective distributions and p is its corresponding probability value. Let M = {m ij } be the cost matrix, where typically m ij = (a ib j ) 2 , then the OT distance can be defined as:

d M (A, B) := min P ∈U (A,B) d i,j=1 p ij m ij (3)
U (A, B) := {P ∈ R d×d + | P 1 d = p A , P T 1 d = p B } (4)
where 1 d is the vector of ones with a dimension of d, and U (A, B) contains all possible joint probabilities of A and B, whose row and column sums to p A and p B respectively. Fig. 3 illustrates an example of OT. The optimal transport plan P , is thus defined as:

P := arg min P ∈U (A,B) d i,j=1 p ij m ij (5)
OT is a special case of the linear programming (LP) problem whose worst case time complexity is O(d 3 logd) [START_REF] Pele | Fast and robust earth mover's distances[END_REF]. Although the acceleration of OT has attracted a large amount of research interest and various algorithms have been proposed, the computational overhead remains high (supercubic). Recently a variant of OT, Sinkhorn distance [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], successfully makes OT several orders of magnitude faster, which greatly enhances the utility of OT. Sinkhorn distance adds an entropic regularization to the classic OT formula and then enforces a simple structure on the optimal transport matrix, as defined by the following equation:

d λ M (A, B) := min P ∈U (A,B)   d i,j=1 p ij m ij + 1 λ d i,j=1 p ij log p ij   ( 6
)
where λ is the regularization coefficient. As λ increases Sinkhorn distance converges to the classic OT distance. When M itself is a metric matrix, Sinkhorn distance satisfies all the three distance axioms, including the symmetry and the triangle inequality. As its name implies, Sinkhorn distance can be computed by Sinkhornâ Ȃ Źs fixed point iteration [START_REF] Franklin | On the scaling of multidimensional matrices[END_REF], [START_REF] Knight | The sinkhorn-knopp algorithm: convergence and applications[END_REF]. To converge to the optimal transport matrix P λ , one only needs to iterate Sinkhornâ Ȃ Źs update a sufficient number of times. The algorithm can be even faster with GPUs because it supports 1-vs-N mode where the distances against multiple targets are computed simultaneously. The 1-vs-N mode turns vector-by-matrix operations into matrixby-matrix operations, which leads to a higher computing density and makes the algorithm more suitable for parallel computing. A flaw of Sinkhorn distance is that sometimes it will be constrained by machine-precision limit when λ increases beyond a problem-dependent value λ max beyond which some elements of e -λM are represented as zeroes.

III. TIME ADAPTIVE OPTIMAL TRANSPORT (TAOT)

The motivation of TAOT is to migrate OT from probability distributions (one-dimensional) to time series data (twodimensional), which means except for the difference in observed values, TAOT also need to capture the difference in time coordinates. An intuitive solution is to consider both observed values and their corresponding time coordinates when calculating the cost between data points in respective time series. In this setting, the one dimensional data point, for example a i , is extended by another time dimension, (a i , t i ), and then a straightforward cost matrix can be defined as

M (i, j) = (a i -b j ) 2 + w * (t i -t j ) 2
, where w is a weight parameter to balance the two parts.

Originally the OT algorithm copes with probability distributions and each input value must have its corresponding probability. Here we assume all observed values of a time series share an equal probability based on the fact that these values come from independent observations. As a result, the input time series can be denoted by:

(A|p A ) = {((a 1 , t 1 )| 1 d ), ((a 2 , t 2 )| 1 d ), ..., ((a d , t d )| 1 d )} (B|p B ) = {((b 1 , t 1 )| 1 d ), ((b 2 , t 2 )| 1 d ), ..., ((b d , t d )| 1 d )} (7
) Since each observed value shares an equal probability, the Sinkhorn iteration [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] can be further simplified. Algorithm 1 summarizes the simplified TAOT algorithm. Basically it contains two phases, where it first prepares the cost matrix and then iterates the simplified Sinkhorn update until the stopping criterion or the maximum iteration number is reached. Note that the time coordinates are normalized into z-score before calculating the cost matrix M . We can get not only the TAOT distance from the algorithm, but also the optimal transport plan, which is significant to the visualization of the results and the execution of other analysises. 24: plan = bsxf un(@times, v , (bsxf un(@times, u, K)))

M (i, j) = (a i -b j ) 2 + w * (t i -t j )
The general idea of TAOT is to integrate temporal dissimilarity between time series to the final distance measure. For Sinkhorn distance based algorithm, the temporal dissimilarity is usually contained in how we construct the cost matrix, and we can observe from Algorithm 1 that any kind of cost matrix can easily fit into this algorithm framework and lead to a new distance measure. In this paper we just employ the most straightforward way to construct our cost matrix. Other definitions of cost matrix is still an open direction for future work.

As for the computational efficiency, given two time series of length N , the time complexity of Euclidean distance is O(N ). The time complexity of DTW-based methods, for example, DTW, derivative DTW, or weighted DTW, is O(N 2 ). The time complexity of TAOT is relatively tricky because it is not as straightforward as the Euclidean distance or most DTW-based methods. TAOT involves a stopping criterion for its iteration and the criterion is calculated based on the current state of the iteration. According to a study on the complexity of multimarginal optimal transport [START_REF] Lin | On the complexity of approximating multimarginal optimal transport[END_REF], the time complexity of TAOT is approximately O(N 2 log N ). Note that TAOT is based on Sinkhorn distance [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], which is a faster version of optimal transport. The naive optimal transport is a computational heavy method, whose cost of computing scales at least in O(N 3 log N ).

IV. EXPERIMENT

A. 1NN ACCURACY REPORT

In this section, we evaluate the classification accuracy of TAOT on 63 datasets of the UCR time series classification archive [START_REF] Chen | The UCR Time Series Classification Archive[END_REF], which has the largest and most widely used collection of public benchmark datasets of time series. These datasets include both real-world time series and synthetic time series from various application domains. The length of time series varies between 60 (Synthetic Control) and 1639 (CinC_ECG_torso), and the number of classes varies between 2 (Gun-Point) and 60 (ShapesAll).

In order to ensure fair comparisons between different similarity measures, we employ the one nearest neighbor classifier (1NN), because the classifier itself does not involve any parameter, and thus the accuracy depends only on the similarity measure. Six most well-established methods, Euclidean distance, DTW, DTW with Sakoe-Chiba constraint, DDTW, WDTW, and FDTW are selected for comparison.

Table 1 summarizes the 1NN classification error rate of TAOT and other competing methods on all datasets. If any parameter is involved, for instance, the width of Sakoe-Chiba band, we choose the optimal one found by cross validation on the training set. From these error rates we can observe that in general TAOT outperforms all the other measures on 28 datasets (Synthetic Control, Swedish Leaf, 50Words, ECG200, Adiac, etc), and achieves leading accuracy together with one or more other measures on another 7 datasets (Wafer, Face(four), Lightning-7, Plane, Coffee, ECGFive-Days, and BirdChicken).

TAOT encounters the aforementioned machine precision limit problem on four datasets (Face(all), Lighting2, Fish, LargeKitchenAppliances) of the original UCR time series classification archive, and therefore there is no accuracy report on these datasets. This is harmless for the classification task since we can find the problem during the early training phase and then change for another alternative method.

Fig. 4(a-b-c-d-e-f) shows the pairwise comparisons between TAOT and the six competing methods, respectively. In Fig. 4, each dot represents a dataset, whose x-coordinate and y-coordinate are respectively the accuracy generated by the competing method and TAOT. In this setting, a dot falling above, on, or below the diagonal indicates that the proposed TAOT outperforms, ties with, or lose to the competing method. The numbers of dots in each different region are labeled on these figures. We can observe that TAOT achieves a better performance in all of the six pairwise competitions. To demonstrate whether TAOT is statistically significant different from other methods, the p-value generated from Friedman test is shown the in lower-right corner of each figure. The pvalue ranges from 0 to 1, and a smaller p-value indicates a more significant difference between methods. Note that all these p-values are smaller than 0.05 and it proves that TAOT performs significantly better than the six competing methods. ) a r is the radius of best Sakoe-Chiba band, measured in percentages of time series length. Optimal values of r are published along with the datasets. b g is a parameter of WDTW, which controls the level of penalization for the points with larger time gap. Optimal values of g are reported by the paper proposing WDTW [START_REF] Jeong | Weighted dynamic time warping for time series classification[END_REF] c Optimal values of λ and w are extracted from cross validations, the details of which will be introduced in Section IV-B.

(a) (b) (c) (d) (e) (f)

B. EXTRACTING RELIABLE PARAMETERS FOR TAOT

TAOT involves two parameters: the regularization coefficient λ and the time weight w. To find a reliable combination of λ and w, we employ grid-search based leave-one-out cross validation on the training set. We choose leave-one-out cross validation over hold-out or k-fold validation because sometimes the training set is too small to further divide. Typically, grid-search is considered to be time-consuming, but for TAOT the search spaces of λ and w are both small and thus the efficiency is acceptable.

For w, if it is too large, for example w > 10, the difference in time coordinates will predominate. Contrarily, if w is too small, for example w < 0.1, the difference in observed values will predominate. Either way we lose the balance. Therefore we set the search space of w to a moderate range, [0.1, 10].

a uniform step size to the and the fraction part respectively, our candidates for w are 0.1, 0.2, ..., 0.9 ..., 10. For λ, recall that as it increases, Sinkhorn distance converges to the classic OT distance. Empirically, 200 is large enough to Sinkhorn distance close to the classic OT distance, and 5 is strict enough for regularizing. So we set the search space of λ to [START_REF] Liao | Clustering of time series data ąła survey[END_REF]200] with a fixed step size of 5.

During grid-search, each pair of λ and w results in an error count, error(λ, w). There can often be multiple pairs of parameters that all achieve the lowest error count, and these optimal parameters on the training set may not perform as well on the testing set. Therefore we develop another strategy to decide the sole and most adaptable pair of parameters. For λ, the most adaptable one should achieve the lowest cumulative error count given every different w, as defined by the following Eq. 8. If there are ties, we choose the smallest λ because it leads to a faster computation.

λ = min arg min λ N n=1 error(λ, w n ) ( 8
)
After λ is fixed, the best w should be the one to have the lowest error count when pairing with the λ , as defined by the following Eq. 9. If there are ties, we choose the median. Fig. 5 shows the real testing error counts of different w candidates that all achieve the lowest training error count. We can observe that the medians of those w candidates achieve the most stable performance in general on the three datasets.

C. THE TEXAS SHARPSHOOTER FALLACY

To further demonstrate the capability of TAOT, we will test whether TAOT is capable of overcoming the Texas Sharpshooter Fallacy. The Texas Sharpshooter Fallacy is a common logic error that happens when comparing methods on multiple datasets. A pervasive scenario is that as long as a method can on some datasets, the author then claim the method is valuable. Because the author think since the method can win on the datasets from some domains, it could be useful in those domains. However, it is not enough to have a method that can be more accurate on some datasets unless you can tell ahead of time that on which datasets it will be more accurate [51]. One way to show whether a method can tell in advance that it can be more accurate on a certain dataset is to test whether the expected accuracy gain over another competing method coincides with the actual accuracy gain. The expected accuracy gain and the actual accuracy gain are defined by Eq. 10 and Eq. 11, respectively. In our setting, the expected accuracy gain is based on the best training accuracy during leave-oneout cross validation, and the actual accuracy gain is based on the best testing accuracy. An expected accuracy gain larger than one indicates that we predict TAOT will perform better, while an actual accuracy gain larger than one indicates that TAOT indeed performs better.

As shown in Fig. 6, Texas sharpshooter plot is a convenient tool to visualize the comparison between the expected accuracy gain and the actual accuracy gain on multiple datasets. Each point represents a dataset to test and each dataset falls into one of the following four possibilities:

• TP(True Positive): In this region we predicted TAOT would increase accuracy, and TAOT did. Obviously, this is the most beneficial situation for TAOT and the majority of points fall into this region. • TN(True Negative): In this region we predicted TAOT would decrease accuracy, and TAOT did. This is not a bad case. Because if we know ahead of time that TAOT will do worse, we can choose another method to avoid the loss of accuracy. • FN(False Negative): In this region we predicted TAOT would decrease accuracy, but the accuracy actually increased. This is also not a bad case. We might miss the opportunity to improve, but we will not do worse. • FP(False Positive): In this region we predicted TAOT would increase accuracy, but the accuracy actually decreased. This is the truly bad case. But not many points fall into this region.

D. THE PAIRING OF MAXIMUM VALUES

One motivation of TAOT is that we want global maximum values to pair with each other if such a pairing is critical for the current problem. Classic OT, which is a special case of TAOT when w = 0, can guarantee the pairing of global maximum values in theory. As for TAOT, if this pairing is decisive, then normally we will get a relatively small w from the training phase, and that will lead to a high possibility of generating a maximum-to-maximum alignment. Fig. 7 gives an example to visually demonstrate that TAOT is more likely to generate the required maximum-to-maximum alignment than the other currently established methods.

E. VISUALIZATION OF ALIGNMENTS GENERATED BY TAOT

The many-to-many alignment generated by TAOT, namely the transport plan, is usually represented by a weight matrix.

As the most established and popular similarity measure of time series, DTW visualizes its output alignment with a warping path in a warping matrix, as shown in the previous Fig. 2. In order to keep consistency with DTW, the visualization of TAOT is still based on a warping matrix, but with a fuzzy warping path or a heat map of the warping path. Fig. 8 illustrates the fuzzy version of a warping path where each

V. CONCLUSION

In this paper, we proposed a new time series similarity measure framework entitled Time Adaptive Optimal Transport (TAOT). It inherits the of several promising properties from the optimal transport (OT) that can greatly alleviate some long-suffered issues with currently established time series similarity measures. To make OT capable of Compared with other currently established methods, TAOT exhibits a better accuracy and it even has the ability to predict whether there will be an accuracy improvement on the majority of experimental datasets. We also introduced a grid-search based parameter extracting strategy and a fuzzy warping path based visualization method for TAOT. TAOT is not just a single distance measure, but also an algorithm framework where other kinds of cost matrix can easily fit into and each creates a new distance measure. We believe this may lead to some potential future work.

FIGURE 1 .

 1 FIGURE 1. Examples of different categories of alignments between time series.

FIGURE 2 .

 2 FIGURE 2. Warping path and cost matrix of DTW.

Algorithm 1

 1 Computation of Time Adaptive Optimal Transport Input: A, B, λ, w, maxIter = 5000, tolerance = 0.005 Output: distance, plan 1: d = |A| 2: t = zscore(linspace(1, d, d)) 3: for i = 1 : |A| do 4: for j = 1 : |B| do 5:

2

 2

8 :

 8 M = M/median(M (:)) 9: K = exp(-λ * M ) 10: curIter = 0 11: u = ones(d, 1)/d 12: while curIter < maxIter do 13: v = 1./(d * K * u) 14: u = 1./(d * K * v) 15: curIter = curIter + 1 16: if mod(curIter, 20) == 1 then 17: criterion = sum(abs(v. * (K * u) -1/d)) 18:if criterion < tolerance then 19:

  end while 23: distance = sum(u. * ((K. * M ) * v))

FIGURE 4 .

 4 FIGURE 4. 1NN accuracy of TAOT versus the six competitors respectively.

VOLUME 4, 2016

 2016 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2020.3016529, IEEE Access Zheng Zhang et al.: Time Adaptive Optimal Transport : A Framework of Time Series Similarity Measure (a) CBF dataset (b) FaceFour dataset (c) Plane dataset

FIGURE 5 .

 5 FIGURE 5. Performance of different candidates of w on testing sets.

  expected gain = training accuracy(T AOT ) training accuracy(competing method) (10) actual gain = testing accuracy(T AOT ) testing accuracy(competing method)

8 VOLUME 4 FIGURE 6 .

 846 FIGURE 6. Texas sharpshooter plot of TAOT versus the other six competing methods respectively.

FIGURE 7 .

 7 FIGURE 7. An example where only TAOT generates the required maximum-to-maximum alignment while the other methods fail.

FIGURE 8 .

 8 FIGURE 8. A fuzzy warping path generated by TAOT.

TABLE 1 .

 1 Comparisons on 1NN Error Rate

	Dataset	Euclidean	DTW	DTW	DDTW	WDTW	FDTW	TAOT
		distance		(r%) a		(g) b		(λ and w) c
	Synthetic Control	0.120	0.007	0.017 (6)	0.433	0.010 (0.30)	0.007	0.003 (195 4)
	Gun-Point	0.087	0.093	0.087 (0)	0.007	0.027 (0.20)	0.020	0.027 (50 0.3)
	CBF	0.148	0.003	0.004 (11)	0.408	0.009 (0.08)	0.001	0.006 (145 1)
	OSU Leaf	0.479	0.409	0.388 (7)	0.120	0.479 (0.60)	0.409	0.450 (160 3)
	Swedish Leaf	0.211	0.208	0.154 (2)	0.115	0.173 (0.03)	0.208	0.107 (175 0.9)
	50Words	0.369	0.310	0.242 (6)	0.308	0.253 (0.10)	0.268	0.233 (150 2)
	Trace	0.240	0.000	0.010 (3)	0.000	0.000 (0.01)	0.000	0.020 (200 0.3)
	Two Patterns	0.090	0.000	0.002 (4)	0.003	0.000 (0.01)	0.000	0.010 (5 6)
	Wafer	0.005	0.020	0.005 (1)	0.022	0.003 (0.30)	0.008	0.003 (70 8)
	Face(four)	0.216	0.170	0.114 (2)	0.375	0.125 (0.10)	0.102	0.102 (40 5)
	Lightning-7	0.425	0.274	0.288 (5)	0.425	0.274 (0.10)	0.301	0.274 (10 0.9)
	ECG200	0.120	0.230	0.120 (0)	0.170	0.130 (0.50)	0.180	0.110 (130 3)
	Adiac	0.389	0.396	0.391 (3)	0.414	0.366 (0.10)	0.414	0.289 (185 0.1)
	Yoga	0.170	0.164	0.155 (2)	0.180	0.153 (0.10)	0.154	0.159 (75 0.2)
	Plane	0.038	0.000	0.000 (6)	0.000	0.000 (0.01)	0.000	0.000 (125 0.5)
	Car	0.267	0.267	0.233 (1)	0.267	0.217 (0.19)	0.367	0.283 (160 0.8)
	Beef	0.333	0.367	0.333 (0)	0.333	0.300 (0.20)	0.367	0.333 (100 6)
	Coffee	0.000	0.000	0.000 (0)	0.071	0.000 (0.01)	0.000	0.000 (85 2)
	OliveOil	0.133	0.167	0.133 (0)	0.133	0.167 (0.01)	0.100	0.133 (200 0.6)
	CinC_ECG_torso	0.103	0.349	0.070 (1)	0.289	0.075 (0.08)	0.317	0.084 (195 10)
	DiatomSizeReduction	0.065	0.033	0.065 (0)	0.065	0.036 (0.10)	0.033	0.020 (10 0.2)
	ECGFiveDays	0.203	0.232	0.203 (0)	0.314	0.138 (0.60)	0.117	0.117 (20 5)
	FacesUCR	0.231	0.095	0.088 (12)	0.157	0.078 (0.03)	0.095	0.067 (30 3)
	ItalyPowerDemand	0.045	0.050	0.045 (0)	0.086	0.043 (0.10)	0.033	0.038 (15 7)
	MedicalImages	0.316	0.263	0.253 (20)	0.337	0.263 (0.08)	0.280	0.296 (40 4)
	MoteStrain	0.121	0.165	0.134 (1)	0.284	0.142 (0.30)	0.165	0.107 (20 1)
	SonyAIBORobotSurface	0.305	0.275	0.305 (0)	0.270	0.255 (0.50)	0.304	0.166 (95 2)
	SonyAIBORobotSurfaceII 0.141	0.169	0.141 (0)	0.142	0.154 (0.05)	0.178	0.134 (105 10)
	Symbols	0.101	0.050	0.062 (8)	0.029	0.049 (0.03)	0.060	0.062 (115 0.8)
	TwoLeadECG	0.253	0.096	0.132 (4)	0.005	0.111 (0.05)	0.112	0.070 (35 0.1)
	Cricket_X	0.423	0.246	0.228 (10)	0.369	0.210 (0.03)	0.269	0.267 (30 3)
	Cricket_Y	0.433	0.256	0.241 (17)	0.441	0.238 (0.01)	0.244	0.277 (145 3)
	Cricket_Z	0.413	0.246	0.254 (5)	0.444	0.246 (0.05)	0.233	0.287 (105 5)
	InsectWingbeatSound	0.438	0.645	0.415 (1)	0.757	0.431 (0.20)	0.591	0.429 (75 10)
	ArrowHead	0.200	0.297	0.200 (0)	0.223	0.183 (0.50)	0.223	0.177 (150 3)
	BeetleFly	0.250	0.300	0.300 (7)	0.250	0.300 (0.01)	0.350	0.100 (10 0.3)
	BirdChicken	0.450	0.250	0.300 (6)	0.150	0.250 (0.01)	0.300	0.150 (45 0.1)
	Ham	0.400	0.533	0.400 (0)	0.524	0.429 (0.20)	0.432	0.371 (40 0.7)
	Herring	0.484	0.469	0.469 (5)	0.500	0.453 (0.01)	0.562	0.297 (60 0.2)
	ProximalPhalanxOAG	0.215	0.195	0.215 (0)	0.180	0.195 (0.01)	0.195	0.185 (15 0.1)
	ProximalPhalanxOC	0.192	0.217	0.210 (1)	0.182	0.213 (0.10)	0.216	0.196 (40 0.7)
	ProximalPhalanxTW	0.293	0.244	0.244 (2)	0.270	0.260 (0.03)	0.288	0.215 (15 0.7)
	ToeSegmentation1	0.320	0.228	0.250 (8)	0.215	0.219 (0.01)	0.276	0.171 (35 0.1)
	ToeSegmentation2	0.192	0.162	0.092 (5)	0.315	0.115 (0.03)	0.154	0.077 (5 0.8)
	DistalPhalanxOAG	0.374	0.230	0.374 (0)	0.240	0.225 (0.40)	0.223	0.185 (5 1)
	DistalPhalanxOC	0.283	0.283	0.275 (1)	0.220	0.237 (0.03)	0.238	0.213 (45 0.4)
	DistalPhalanxTW	0.367	0.410	0.367 (0)	0.273	0.268 (0.30)	0.278	0.245 (5 0.5)
	Earthquakes	0.288	0.281	0.273 (6)	0.276	0.292 (0.03)	0.276	0.174 (5 7)
	MiddlePhalanxOAG	0.481	0.500	0.481 (0)	0.255	0.260 (0.05)	0.283	0.238 (60 0.2)
	MiddlePhalanxOC	0.234	0.302	0.234 (0)	0.278	0.292 (0.30)	0.257	0.227 (20 0.5)
	MiddlePhalanxTW	0.487	0.494	0.494 (3)	0.444	0.414 (0.08)	0.416	0.381 (70 0.4)
	ShapeletSim	0.461	0.350	0.300 (3)	0.461	0.244 (0.03)	0.122	0.439 (60 2)
	Wine	0.389	0.426	0.389 (0)	0.481	0.426 (0.10)	0.370	0.296 (55 9

TABLE 2 .

 2 Comparisons on 1NN Error Rate (Continued)

	Dataset	Euclidean	DTW	DTW	DDTW	WDTW	FDTW	TAOT
		distance		(r%) a		(g) b		(λ and w) c
	WordsSynonyms	0.382	0.351	0.262 (9)	0.320	0.249 (0.05)	0.343	0.306 (55 3)
	Computers	0.424	0.300	0.380 (12)	0.332	0.416 (0.01)	0.348	0.340 (15 0.6)
	Meat	0.067	0.067	0.067 (0)	0.333	0.067 (0.10)	0.067	0.017 (95 0.9)
	RefrigerationDevices	0.605	0.536	0.560 (8)	0.592	0.592 (0.03)	0.536	0.469 (15 0.2)
	ScreenType	0.640	0.603	0.589 (17)	0.573	0.589 (0.01)	0.603	0.565 (5 0.2)
	ShapesAll	0.248	0.232	0.198 (4)	0.165	0.192 (0.05)	0.232	0.193 (165 0.8)
	SmallKitchenAppliances	0.659	0.357	0.328 (15)	0.347	0.347 (0.01)	0.333	0.336 (35 4)
	Strawberry	0.054	0.059	0.054 (0)	0.059	0.062 (0.40)	0.054	0.060 (95 0.2)
	Worms	0.545	0.416	0.468 (9)	0.497	0.552 (0.03)	0.514	0.575 (35 1)
	WormsTwoClass	0.390	0.377	0.416 (7)	0.298	0.376 (0.05)	0.354	0.359 (25 8)

Zheng Zhang et al.: Time Adaptive Optimal Transport : A Framework of Time Series Similarity Measure

 VOLUME 4, 2016

 VOLUME 4, 2016

 VOLUME 4, 2016

 VOLUME 4, 2016

This work was supported by the National Natural Science Foundation of China (Grant No.41701399), and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No.XDA19080301).

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2020.3016529, IEEE Access

Currently, she is serving as a Team Leader and leading a project of multispectral imagery radiometric and geometric correction for large volume of image data at global scale for higher resolution global land-cover mapping. She has significant experience in managing and designing software systems for satellite image processing and applications. Her research interests include the use of mathematical theories to development algorithms related to image processing and analysis.

THOMAS CORPETTI received an engineer degree in electrical engineering and a master degree in computer vision, in 1999. He obtained his Ph.D. and habilitation degree in 2002 and 2011, in the field of Computer Vision and Applied Mathematics with University Rennes I, France. After Ph.D. he spent one year as an assistant professor in Rennes and one year as a post-doc with IRSTEA (environment institute). He then obtained a permanent researcher position at CNRS in 2004 (French National Institute for Scientific Research) on the analysis of remote sensing image sequences for environmental applications. From 2009 to 2012, he has been with LIAMA, a sino-french lab in computer sciences, automatics and applied mathematics at Beijing, China, where he headed the TIPE group (Turbulence, Images, Physics and Environment).

He is now with OSUR (Observatory for Universe Sciences of Rennes, France) and LETG (Littoral, Environnement, Télédétection, Géomatique) UMR 6554 as a director of research, CNRS. His main research interests are concerned with the definition of computer vision tools for the analysis of remote sensing data (low and high resolution) for environmental applications.