\section{Introduction}\label{sec:Intro}

With constant increase of computational power data-driven methods have gained a significant popularity in the past few decades with applications spanning all fields of science \cite{Montans19}. The study of {\it intrinsic localized modes} (ILM) in crystal lattice models, also known as {\it discrete breathers} (DB), is an active field of research in mathematics and physics \cite{MacKay94,Flach08,ProcMica14,ram20} which relies heavily on obtained data from numerical simulations. Thus, emerging data-driven methods provide great opportunity to explore new techniques for understanding nonlinear wave phenomena in crystal lattice models.  

Such localized wave solutions are of particular interest in material science where transport of charge in silicates by moving nonlinear localized modes is experimentally confirmed, the phenomenon called {\it hyperconductivity} \cite{ram20,Russell19}. In addition, localized modes provide mechanisms for the formation of long decorated dark lines in muscovite mica \cite{ProcMica14} and can have a functional role in proteins \cite{Senet15,Piazza08}. Strong evidence of localized phenomena in crystal lattice models is provided by numerical experiments \cite{Chris11,Marin01,Bajars15,Shepelev20,Chetverikov16,Bajars21}, to mention but a few.

There are already well established theoretical and numerical understandings of the existence of one-dimensional stationary DB \cite{MacKay94,Flach08,Aubry06}, i.e., spatially localized time-periodic excitations. The same may not be said about propagating DB \cite{Flach99} and, in particular, in higher dimensions. There are still open theoretical questions regarding the existence and properties of propagating DB in general nonintegrable lattices. Mobile discrete breathers are not expected to be long-lived, since the lattice models considered are likely to be nonintegrable and their lifespan is subject to interactions with lattice defects and the phonon background, which may be viewed as thermal noise. Thus, an obvious challenge is to study the existence and interactions of long-lived propagating discrete breather solutions \cite{Bajars15,Bajars21} to understand their role in the transport of energy in sputtering experiments \cite{Russell07}. Such demanding computations require good data-based numerical tools to detect and track breather solutions in space and time which serves as a motivation for our work. 

In this work we develop classification algorithms to differentiate between nonlinear localized waves and nonlocalized linear waves in numerical simulations of one-dimensional crystal lattice model. Obtained classifiers rely on locally sampled data in opposed to nonlocal time series analysis methods, such as discrete Fourier and wavelet transforms. Such classifiers can be efficiently trained with different given labeled datasets which are not only limited to numerical simulation data. Trained classifiers then can be used to detect localization regions, e.g., in numerical simulations, which is the first step towards fully automated tool for quantitative data-based analysis of complex numerical experiments. Importantly, our analysis and methodology, which follows, extends, in general, to any one-dimensional crystal lattice models which support intrinsic localized mode solutions.    

The paper is organized in the following way. In Section \ref{sec:MathModel} mathematical model of one-dimensional crystal lattice is formulated. Spectral properties of discrete lattice waves are discussed in Section \ref{sec:Spectrum}. Dimensionality reduction methods and classification algorithms are presented in Section \ref{sec:Classification}. In Section \ref{sec:Applications} application of the classification algorithms is demonstrated for detecting localized modes in numerical simulations. Discussion and conclusions are provided in Section \ref{sec:Conclusions}.    

\section{Mathematical model}\label{sec:MathModel}

In this work we are concerned with one-dimensional crystal lattice models of classical molecular dynamics exhibiting ILM solutions. In particular, we consider Hamiltonian dynamics of $N$ particles arising from the dimensionless Hamiltonian
\begin{equation}\label{eq:Hamilt}
H = K + U + V 
  = \sum_{n=1}^{N} \left( \frac{1}{2} \dot{q}_{n}^2 +
  U(q_{n}) + V(|q_{n+1}-q_{n}|) \right),
\end{equation}
where $K$ is the kinetic energy, $U$ is the on-site potential energy, $V$ is the interaction potential of particles, $q_n\in\mathbb{R}$ is the position of the $n^{th}$ particle, $\dot{q}_{n}$ is its time derivative and $|\cdot|$ is the Euclidean distance. In the dimensionless form \eqref{eq:Hamilt} particle equilibrium distance $\Delta=q_{n+1}-q_n$ as well as all masses are equal to one. 

For example, such model \eqref{eq:Hamilt} arises as one-dimensional model of the K-K layer of layered silicate muscovite mica \cite{Chris11} or simplified model of two-dimensional hexagonal cristal lattice \cite{Bajars15,Bajars21}. In such case, the on-site potential $U$ models forces from atoms above and below the K-K sheet. Thus, the on-site potential is modeled as smoothed periodic function \cite{Chris11}:
\begin{equation}\label{eq:OnSite}
U(q_n) = 1-\cos(2 \pi q_n).
\end{equation}
Note that more realistic potentials for muscovite mica can also be considered \cite{Juan19}.

There are multiple options for modeling the interaction potential $V$ \cite{Allen89}. In this work we consider well known scaled Lennard-Jones potential in the following form:
\begin{equation}\label{eq:LJ}
V(r_n) = \epsilon \left( \left( \frac{1}{r_n} \right)^{12} - 2 \left(
    \frac{1}{r_n} \right)^{6} \right), \quad r_n = |q_{n+1}-q_n|,
\end{equation}
where the parameter $\epsilon>0$ describes the ratio of the atomic well depths associated to the interaction and the on-site potential. That is, large values of $\epsilon$ lead to a model describing one-dimensional Lennard-Jones fluid, while small values of $\epsilon$ lead to system of decoupled nonlinear oscillators. Further simplifications in \eqref{eq:Hamilt} have been made to only consider close neighbor interactions with periodic boundary conditions, i.e., ${q}_{N+1} = {q}_{1}$.

From the Hamiltonian \eqref{eq:Hamilt} obtained Hamiltonian equations:
\begin{align}
\dot{q}_n & = p_n, \label{eq:q}\\
\dot{p}_n & = -U'(q_n) + V'(|q_{n+1}-q_{n}|) - V'(|q_{n}-q_{n-1}|), \label{eq:p}
\end{align}
where $p_n$ is momentum, are solved numerically with the second order time reversible symplectic Verlet method \cite{Hairer}. In the following, without loss of generality, all analysis and numerical simulations are performed with time step $\tau=0.01$ and $\epsilon=0.05$. 

To excite DB solutions we consider particles in their dynamical equilibrium states, i.e., $q^0_n=n-1$ and $p^0_n=0$, while exciting four neighboring particle momenta with the pattern
\begin{equation}\label{eq:patern}
\B{p}_0 = \gamma \, (-1, 2, -2, 1)^T, \quad \gamma > 0.
\end{equation} 
Such excitation, due to nonlinearity, generates in our model stationary discrete breather. In the current study we only consider numerical data from stationary DB solutions while the extension to differentiate between stationary and moving breather solutions is left for future work.

In general, larger values of $\gamma$ generate localized waves with larger particle displacements. While such excitation \eqref{eq:patern} does not generate exact wave solutions the amount of generated linear phonon waves has been observed to be minimal. Not all values of $\gamma$ can be considered. Too small values will generate phonons, while too large values will force particles to move out from their potential wells or in addition to the stationary breather will also generate moving breathers. In addition, appropriate values of $\gamma$ depend on $\epsilon$. With spectral analysis of Section \ref{sec:Spectrum} and $\epsilon=0.05$ we find $\gamma$ lower bound $0.15$ for the initial momenta pattern \eqref{eq:patern} to generate stationary DB. On the other hand, to generate arbitrary phonon waves we can choose small random momenta values $p_n$ as initial conditions, as confirmed by the spectral analysis in Section \ref{sec:Spectrum}. 

Both described initial conditions will be used to generate datasets in Section \ref{sec:Classification} to obtain classification and detection algorithms which are capable of distinguishing between localized and phonon waves. To see discrepancies between both types of waves in Figure \ref{fig:sol} we illustrate particle displacement function $u_n = q_n - q_n^0$ from two numerical simulations with $N=24$ at final computational time $T_{end}=5$. In Figure \ref{fig:Bsol} we show displacement function of a localized discrete breather solution computed with $\gamma=0.3$, while in Figure \ref{fig:Psol} we demonstrate particle displacements for phonon waves computed with random initial momentum values. It is easy to see larger displacement values for the breather solution and that it is spatially localized compared to the approximately linear waves. In addition, in Figure \ref{fig:Bsol} we have added markers to eight particle displacement values on which the breather is localized. Exactly eight data point values will be used in the classification algorithms of Section \ref{sec:Classification}.

\begin{figure}[t]
\centering 
\subfigure[]{\label{fig:Bsol}
\includegraphics[trim=1cm 0.1cm 1cm 1cm,clip=true,width=0.48\textwidth]{Fig1a.jpg}}
\subfigure[]{\label{fig:Psol}
\includegraphics[trim=1cm 0.1cm 1cm 1cm,clip=true,width=0.48\textwidth]{Fig1b.jpg}}
\caption{Particle displacements from the numerical simulations at $T_{end}=5$ with $\epsilon=0.05$ and $N=24$. (a) stationary discrete breather solution initiated with $\gamma=0.3$ and markers indicating eight particles on which the breather is localized. (b) phonon wave solution which does not exhibit localization.}\label{fig:sol}
\end{figure}

\section{Spectral properties of lattice waves} \label{sec:Spectrum}

It is well known that linear phonon waves in crystal lattice models can be characterized by their dispersion relation \cite{Flach08}. For the model \eqref{eq:q}--\eqref{eq:p} the dispersion relation is
\begin{equation}\label{eq:dispersion}
\omega^2 = \omega_0^2 + 4 c^2 \sin^2 \left( \frac{k}{2} \right),
\end{equation}
where $\omega$ and $k$ are the phonon frequency and wave number, respectively, $\omega_0^2=\frac{\partial^2 U}{\partial q_n^2}(q_n^0)=4\pi^2$ is the oscillation frequency of the isolated oscillators when $\epsilon=0$ and $ c^2 = \frac{\partial^2 V}{\partial q_n^2}(q_n^0)= 72 \epsilon $ is the speed of sound.

In Figure \ref{fig:SP} we plot normalized two-dimensional spectrum of numerical simulation data in space and time, i.e., contours of 2D discrete Fourier transform. We consider a lattice of $N=128$ particles and integrate in time until $T_{end}=400$. The solid line indicates the dispersion relation \eqref{eq:dispersion}. 

Excited stationary breather solution by \eqref{eq:patern} are optical breathers with particle frequencies above the phonon linear spectrum as can be seen in Figure \ref{fig:SPb}. Only small amount of energy is localized on the dispersion curve indicating minimal presence of phonons in numerical simulation. 

In Figure \ref{fig:SPp} we considered random initial conditions for the momentum variable $p_n$ drawn from the normal distribution with zero mean and variance $\sigma^2=0.001$. Such initial condition leads to full spectrum of phonons which coincides with the dispersion curve \eqref{eq:dispersion}. Thus, this provides representative data of phonons for data analysis. 

\begin{figure}[t]
\centering 
\subfigure[]{\label{fig:SPb}
\includegraphics[trim=0.5cm 0.5cm 0.8cm 1cm,clip=true,width=0.48\textwidth]{Fig2a.jpg}}
\subfigure[]{\label{fig:SPp}
\includegraphics[trim=0.5cm 0.5cm 0.8cm 1cm,clip=true,width=0.48\textwidth]{Fig2b.jpg}}
\caption{Contours of normalized two-dimensional discrete Fourier transform in space and time of numerical simulation with $N=128$ and $T_{end}=400$. Dispersion relation \eqref{eq:dispersion} is indicated by the solid line. (a) stationary discrete breather spectrum computed with the initial pattern \eqref{eq:patern} and $\gamma=0.3$. (b) phonon wave spectrum computed with random initial condition and coinciding with the dispersion curve \eqref{eq:dispersion}.}\label{fig:SP}
\end{figure}

Spectral analysis of this section not only provides insights into the lattice wave properties but also allows to justify the methods for obtaining a so called labeled training and testing data for classification problem presented in the following section.

\section{Classification of lattice waves} \label{sec:Classification}

In this section we describe classification algorithms of crystal lattice waves based on locally sampled data. Machine learning techniques \cite{ML} are considered for data dimensionality reduction and building an efficient classifier. All calculations are performed in {\it Python} with open source {\it scikit-learn}\footnote{\url{https://scikit-learn.org/stable/}} library.

We differentiate localized discrete breather solutions, see Figure \ref{fig:Bsol}, from nonlocalized linear phonons, see Figure \ref{fig:Psol}. Training (including validation) and testing data is obtained from total of $N_{sim}=1000$ numerical simulations with $0.75$ localized and nonlocalized wave data proportion. Local data of particle displacements, momenta and energies are sampled from eight neighboring particles indicated by red dots in Figure \ref{fig:Bsol}. We consider three different datasets for classification: $X_{u}$, $X_{u,p}$ and $X_{u,p,E}$, where the dataset $X_{u}$ only contains particle displacement $u_n$ values, the dataset $X_{u,p}$ contains particle displacement $u_n$ and momenta $p_n$ values, while the dataset $X_{u,p,E}$ in addition to displacement and momenta values contains also particle energy density data defined by
\begin{equation}\label{eq:En}
E_n := \frac{1}{2} {p}_{n}^2 + U(q_{n}) + \frac{1}{2}\left( V(|q_{n}-q_{n-1}|) + V(|q_{n+1}-q_{n}|) \right).
\end{equation}

We excite discrete breather solutions with momenta pattern \eqref{eq:patern} where parameter $\gamma$ values are randomly chosen from the uniform distribution $\mathcal{U}(0.15,1.2)$. Linear phonon waves are generated with random initial particle displacement and momenta values from the uniform distribution $\mathcal{U}(-0.01,0.01)$.    

\subsection{Dimensionality reduction}

To increase the efficiency of the classification algorithms we consider two dimensionality reduction techniques, i.e., {\it Principal Component Analysis} (PCA) and {\it Locally Linear Embedding} (LLE) \cite{ML}. PCA is a data projection algorithm onto so called principal components. In our experiments we find that with projection down to only two principal components we are able to preserve more than $95\%$ of datasets's variance. On the other hand, LLE is a dimensionality reduction technique that computes low dimensional, neighborhood preserving embeddings of high dimensional data. Performing LLE we considered twenty neighbors for each data point.

In Figure \ref{fig:dimred} we illustrate dimension reduction of all three datasets $X_{u}$, $X_{u,p}$ and $X_{u,p,E}$ applying both dimensionality reduction techniques PCA and LLE, respectively. Colorbars in Figure \ref{fig:dimred} indicate $\gamma$ values of the pattern \eqref{eq:patern} used to obtain particular DB solution data. Apart from the PCA projection of the dataset $X_{u,p}$, see Figure \ref{fig:PCAWP}, the $\gamma$ values can be traced along the continuous curve in two dimensions. Interestingly, in Figure \ref{fig:LLEX} the curve is a straight line parallel to $y$-axis.  

In Figures \ref{fig:PCAW}, \ref{fig:PCAWP} and \ref{fig:PCAX} we projected datasets down onto two principal components $(PC1,PC2)$ while in Figures \ref{fig:LLEW}, \ref{fig:LLEWP} and \ref{fig:LLEX} we show local data linear embedding in two-dimensional $(LLE1,LLE2)$ coordinates. Notice that linear wave data is hardly separated from nonlinear localized wave data in Figures \ref{fig:PCAW}, \ref{fig:LLEW} and \ref{fig:PCAWP}. For the dataset $X_{u,p}$ Figure \ref{fig:LLEWP} suggests that LLE is superior over PCA for the classification problem. Locally linear embedding of the dataset $X_{u,p,E}$ is able to produce the largest separation between two types of wave data. Thus, information of particle displacements and momenta alone may not be sufficient to build an optimal classifier for detecting ILM, e.g., in numerical simulations, see Section \ref{sec:Applications}.     

\begin{figure}[!ht]
\centering 
\subfigure[]{\label{fig:PCAW}
\includegraphics[width=0.47\textwidth]{Fig3a.jpg}}
\subfigure[]{\label{fig:LLEW}
\includegraphics[width=0.47\textwidth]{Fig3b.jpg}}
\subfigure[]{\label{fig:PCAWP}
\includegraphics[width=0.47\textwidth]{Fig3c.jpg}}
\subfigure[]{\label{fig:LLEWP}
\includegraphics[width=0.47\textwidth]{Fig3d.jpg}}
\subfigure[]{\label{fig:PCAX}
\includegraphics[width=0.47\textwidth]{Fig3e.jpg}}
\subfigure[]{\label{fig:LLEX}
\includegraphics[width=0.47\textwidth]{Fig3f.jpg}}
\caption{Dimensionality reduction with PCA and LLE of three datasets $X_{u}$, $X_{u,p}$ and $X_{u,p,E}$ and decision boundaries of SVC. (a) PCA of $X_{u}$ with nonlinear SVC. (b) LLE of $X_{u}$ with nonlinear SVC. (c) PCA of $X_{u,p}$ with nonlinear SVC. (d) LLE of $X_{u,p}$ with nonlinear SVC. (e) PCA of $X_{u,p,E}$ with linear SVC. (f) LLE of $X_{u,p,E}$ with linear SVC.}\label{fig:dimred}
\end{figure}

\subsection{Classification algorithms}

Linear and nonlinear kernel {\it Support Vector Machine} \cite{ML} classifiers (SVC) are trained for lattice wave data classification. SVC are supervised learning algorithms that construct a decision boundary with the largest distance to the nearest training-data point of each data class. The Gaussian {\it Radial Basis Function} (RBF) kernel is used in nonlinear SVC which showed improved performance over other nonlinear kernels, such as polynomial and hyperbolic tangent. Soft margin SVC are considered where optimal hyperparameter values are found with {\it grid search} technique.

To evaluate the performance of linear and nonlinear classifiers for each dataset $\mathrm{K}$-fold cross-validation was considered with $0.7$ training and testing data proportion. In our experiments we considered $\mathrm{K}=5$. The {\it precision} and {\it recall} scores were recorded on cross-validation and testing datasets. The precision is a metric of accuracy of the positive (localized wave) predictions while the recall is a metric of positive instances that are correctly detected by the classifier.

Averaged precision and recall values over $100$ random {\it train-test-splits} of datasets $X_{u}$, $X_{u,p}$ and $X_{u,p,E}$ are listed in Table \ref{tab:Precision_recall}. For each train-test-split we compute precision and recall scores for the testing data and validation dataset obtained from the $\mathrm{K}$-fold of training data. Notice that in most cases recall values are smaller than precision values which is attributed to the fact that more localized wave data is classified as linear opposed to linear waves being classified as localized. In addition to mean values we also computed (not shown) variance values which in all cases, if not exactly zero, did not exceed $10^{-2}$ value.

Examining Table \ref{tab:Precision_recall} we observe exceptionally good results except for the PCA of the dataset $X_{u,p}$ when linear SVC is used. That is also evident in Figure \ref{fig:PCAWP}, where it is easy to see that nonlinear waves can not be separated from linear waves by a straight line. In that case the use of nonlinear SVC is required. Very good results are demonstrated by PCA and LLE dimensionality reduction of the dataset $X_{u,p,E}$, either with linear or nonlinear SVC. In this case for simplicity we advocate the use of linear SVC. Precision and recall values equal to one state that none of the data samples were incorrectly classified by SVC. As follows we consider linear SVC for the dataset $X_{u,p,E}$ and nonlinear classifier for the datasets $X_{u}$ and $X_{u,p}$, which is illustrated in Figure \ref{fig:dimred} where the prediction boundaries and contours of the decision function are drawn. 

In summary, Table \ref{tab:Precision_recall} shows that the classifier is not sensitive to random train-test-split of the datasets and that for the dataset $X_{u,p,E}$ SVC is essentially a perfect classifier for both dimensionality reduction techniques. To further differentiate the classifiers we consider their classification performance on imperfect data, see the following section.           

\begin{table}[t]
	\begin{center}
		\begin{tabular}{|c||c|c||c|c||c|}
			\hline
			\multirow{2}{*}{} & \multicolumn{2}{c||}{{Linear kernel}} &
			\multicolumn{2}{c||}{{RBF kernel}} & \\
			\cline{2-5}
			& {PCA} & {LLE} & {PCA} & {LLE} & \\			
			\specialrule{.1em}{.05em}{.05em} 
			{Precision} (validation set) & 0.9985 & 0.9623 & 1 & 0.9953
			& \multirow{4}{*}{$X_u$} \\
			{Recall} (validation set) & 0.9533 & 0.9478 & 0.9998 & 0.9354 & \\
			\cline{1-5}
			{Precision} (testing set) & 0.9984 & 0.9658 & 1 & 0.9965 & \\
			{Recall} (testing set) & 0.9516 & 0.9487 & 0.9992 & 0.9370 & \\
			\specialrule{.1em}{.05em}{.05em} 
			{Precision} (validation set) & 0.7947 & 0.9927 & 1 & 1
			& \multirow{4}{*}{$X_{u,p}$} \\
			{Recall} (validation set) & 
0.9544 & 0.9993 & 0.9754 & 0.9990 & \\
			\cline{1-5}
			{Precision} (testing set) & 0.7978 & 0.9928 & 1 & 0.9999 & \\
			{Recall} (testing set) & 
0.9519 & 0.9999 & 0.9758 & 0.9991 & \\
						
			\specialrule{.1em}{.05em}{.05em} 
			{Precision} (validation set) & 1 & 1 & 1 & 1 & \multirow{4}{*}{$X_{u,p,E}$} \\
			{Recall} (validation set) & 1 & 1 & 1 & 1 & \\
			\cline{1-5}
			{Precision} (testing set) & 1 & 1 & 1 & 1 & \\
			{Recall} (testing set) & 1 & 1 & 1 & 1 & \\
			\hline
		\end{tabular}
	\end{center}
	\caption{Precision and recall metrics averaged over $100$ random train-test-splits of different classification algorithms, dimensionality reduction techniques and datasets.}
	\label{tab:Precision_recall}
\end{table}

\subsection{Classification of imperfect data}

So far we have considered datasets of perfectly sampled data of linear and localized nonlinear waves. With application in mind, i.e., detecting ILM in numerical simulations, we consider additional datasets $X_{u,p,E}^{n_{off}}$ of localized discrete breather data from $N_{sim}=1000$ simulations with randomly chosen $\gamma$ values from the uniform distribution $\mathcal{U}(0.15,1.2)$. We define index $n_{off}$, where $n_{off}=0,1,\dots,8$, which indicates the shift in particles from which the data is sampled, i.e., the data is sampled from particles $q_{m+n_{off}}$, where particle $q_m$ is the first particle in Figure \ref{fig:Bsol} identified with red dot. Thus, the sampled datasets $X_{u,p,E}^{n_{off}}$ consist of partially linear and nonlinear wave data. 

\begin{figure}[t]
\centering 
\includegraphics[width=0.6\textwidth]{Fig4.jpg}
\caption{Percentages of classified imperfect localized wave data of the datasets $X_{u,p,E}^{n_{off}}$ by linear SVC with dimensionality reduction methods PCA and LLE.}
\label{fig:dataoff}
\end{figure}      

We classify samples from the datasets $X_{u,p,E}^{n_{off}}$ using linear SVC and both dimensionality reduction methods PCA and LLE. In Figure \ref{fig:dataoff} we plot the percentage of wave data being classified as localized wave or linear wave depending on the shift index $n_{off}$. When $n_{off}=0$, the data consists of purely nonlinear wave data and we can see $100\%$ identification of it. On the other hand, when $n_{off}=8$ the dataset does not contain information of nonlinear waves and we see $100\%$ identification of linear wave data. Both classifiers demonstrate gradual decrease in localized wave identification as $n_{off}$ value increases. Notice that the curve of LLE drops faster to zero compared to PCA. This has an impact on detection region size, i.e., SVC with PCA produces wider detection regions of localization as indicated by numerical simulations of Section \ref{sec:Applications}.
 
\section{Applications}\label{sec:Applications}

In this section we apply classification algorithms developed above to detect localized DB of the one-dimensional crystal lattice model \eqref{eq:Hamilt}. We consider sliding window object detection method by sliding a window of eight particle size at each time step over the computational data. To detect region of a localized wave we define a discrete density function $\rho_n$, $n=1,\dots,N$, which takes a positive value if data at particular eight particles is classified as being of localized wave, otherwise the density value is equal to zero. In the following Figures \ref{fig:ApplA} and \ref{fig:ApplBC} we identify all particles by red dots where the localization is detected.

For the first example we consider a lattice of $N=64$ particles and excite two stationary breathers using momenta pattern \eqref{eq:patern} with $\gamma=0.55$ and $\gamma=0.35$ at particles with indices $(16,17,18,19)$ and $(44,45,46,47)$, respectively. We solve the Hamiltonian system \eqref{eq:q}--\eqref{eq:p} until time $T_{end}=1$ and detect the regions of localization, see Figure \ref{fig:ApplA}.

\begin{figure}[t]
\centering 
\subfigure[]{\label{fig:ApplA_1}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig5a.jpg}}
\subfigure[]{\label{fig:ApplA_2}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig5b.jpg}}
\subfigure[]{\label{fig:ApplA_3}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig5c.jpg}}
\subfigure[]{\label{fig:ApplA_4}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig5d.jpg}}
\caption{Detection of numerically simulated two stationary breathers excited by momenta pattern \eqref{eq:patern} with $\gamma=0.55$ and $\gamma=0.35$ on the left and right, respectively. Detected localization region particles are identified by red dots. (a) linear SVC with PCA trained on the dataset $X_{u,p,E}$. (b) linear SVC with LLE trained on the dataset $X_{u,p,E}$. (c) nonlinear SVC with PCA trained on the dataset $X_{u}$. (d) nonlinear SVC with LLE trained on the dataset $X_{u}$.}\label{fig:ApplA}
\end{figure}

For this experiment we used four classifiers. In Figure \ref{fig:ApplA_1} we demonstrate results of linear SVC classifier with PCA, while in Figure \ref{fig:ApplA_2} we consider linear SVC classifier with LLE. Both classifiers are trained on the dataset $X_{u,p,E}$ and are able to capture regions of localization very well. As previously suggested, see Figure \ref{fig:dataoff}, on average SVC classifier with PCA produces wider detection regions, as can be seen from number of colored particles.

In comparison, nonlinear classifiers trained on the dataset $X_{u}$ containing only the information of particle displacements, see Figures \ref{fig:ApplA_3} and \ref{fig:ApplA_4}, are not always able to capture small displacement ILM. Similar issues (not shown) were observed applying nonlinear classifiers trained on the dataset $X_{u,p}$. Despite these shortcomings the classifiers can still be applied if additional information of the energy \eqref{eq:En} is unknown.

To further demonstrate methods' robustness and future research directions we consider two additional examples. In the first example we consider a stationary discrete breather in noisy background by initially randomly exciting all particles and their momenta with addition to the pattern \eqref{eq:patern} with $\gamma=0.5$. In this example we performed calculations until time $T_{end}=10$ and used linear SVC classifier with LLE trained on the dataset $X_{u,p,E}$, see Figure \ref{fig:ApplC}. Figure \ref{fig:ApplC} illustrates very well that the method was able to detect the regions of localization despite the presence of large amount of phonon waves. These results suggest that the classification algorithms can be applied to thermostatic dynamics, i.e., during experiments at fixed temperature.

\begin{figure}[t]
\centering 
\subfigure[]{\label{fig:ApplC}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig6a.jpg}}
\subfigure[]{\label{fig:ApplB}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig6b.jpg}}
\caption{Detection of localization regions in numerical simulations with linear SVC and LLE trained on the dataset $X_{u,p,E}$. Particles of detected regions are identified by red dots. (a) simulation of stationary discrete breather ($\gamma=0.5$) in noisy background. (b) simulation of two mobile breather collision initiated by momenta pattern \eqref{eq:pattern2} with $\gamma=0.45$ (left) and $\gamma=-0.65$ (right).}\label{fig:ApplBC}
\end{figure}

In the second example we consider a simulation of two mobile breather collision initiated by three neighboring particle momenta with pattern
\begin{equation}\label{eq:pattern2}
\B{p}_0 = \gamma \, (-1, 2, -1)^T.
\end{equation}
The mobile breather propagating to right was excited with $\gamma=0.45$ while the mobile breather propagating to left was excited with $\gamma=-0.65$. Integration in time was performed until $T_{end}=100$. Around time $t=50$ both breathers collide and pass through or reflect from each other and continue to propagate, see Figure \ref{fig:ApplB}. Although the classifier was trained on data obtained from stationary discrete breathers, we are still able to detect the localization regions before, during and after the collision. As in Figure \ref{fig:ApplC} we used linear SVC with LLE trained on the dataset $X_{u,p,E}$. These satisfactory results show that the methodology can be used beyond just detection of stationary DB.

\section{Discussion and conclusions}\label{sec:Conclusions}

The results of Section \ref{sec:Applications} demonstrate that classification algorithms, such as SVC, either linear or nonlinear, can be used to detect localization regions in numerical simulations of ILM. The strength of the developed approach lies in the fact that localization regions can be detected from locally sampled data. In future research we will look for the method which makes predictions with the most optimal width of the localization region. For example, in this work we only considered data sampled from eight neighboring particles, other options could also be investigated.

Numerical examples demonstrated that our classifiers are also able to detect mobile DB but lack the ability to differentiate between these two types of localized solutions. To address this we will explore multilabel and multioutput classification algorithms in future work. In addition, we plan to extend our work to higher dimensional crystal lattice models and include kink solutions.

For efficiency and visualization purposes with high success we applied two dimensionality reduction methods before training SVC. The best results were obtained with LLE on the dataset $X_{u,p,E}$, such that we could use linear SVC and find more optimal localization regions in numerical experiments. In future work we plan to explore, in particular, more other manifold learning techniques.  

To demonstrate versatility of our approach we considered three different datasets for training. We found that the best results are obtained if in addition to particle displacement values particle momenta and energy density values are also added to the dataset. Nevertheless, classifiers can be built on all three datasets and overall good performance can be achieved.


\section*{Acknowledgements}
J.~Baj\={a}rs acknowledges support from PostDocLatvia grant No.1.1.1.2/VIAA/4/20/617.

\bibliographystyle{unsrt}  
