Annotation and Detection of Emotion in Text-based Dialogue Systems with CNN

Knowledge of users' emotion states helps improve human-computer interaction. In this work, we presented EmoNet, an emotion detector of Chinese daily dialogues based on deep convolutional neural networks. In order to maintain the original linguistic features, such as the order, commonly used methods like segmentation and keywords extraction were not adopted, instead we increased the depth of CNN and tried to let CNN learn inner linguistic relationships. Our main contribution is that we presented a new model and a new pipeline which can be used in multi-language environment to solve sentimental problems. Experimental results shows EmoNet has a great capacity in learning the emotion of dialogues and achieves a better result than other state of art detectors do.

Introduction
Emotion detection is potentially important in voice recognition systems, which provides improvements to further humancomputer interactions.
Emotion detection can be achieved by various approaches, including text, gestures, facial expression, electrocardiographs, etc [1][2][3][4]. Much work based on these features has been introduced over the last decade [5]. However, text is the most common way for communication in the Internet. Textbased emotion detection has multiple applications in humancomputer interaction, especially in search engines and chatting robots. With the knowledge of users' emotion, search engines could give more specific results, while chatting robots could reply in a more natural way.
Methods in text-based emotion detection can be divided into three categories: keyword-based detection, learning-based detection, and hybrid detection [5]. CNNs (convolutional neural networks) have shown their great capacity of learning the nature of images and speeches over the past few years. With a deep configuration and a large training dataset, keyword extraction might be not necessary anymore, which will significantly simplify the classification process.
In this work, we presented and tested a Chinese daily dialogues emotion detector based on deep convolutional neural network without segmentation or keyword extraction. This paper is divided into 7 sections. Section 2 gives an introduction of the pre-processing step. Section 3 gives the CNN configuration of EmoNet. Section 4 gives details of the training process. In section 5, we tested different configurations and hyperparameters and gave evaluations based on top-1 accuracy and comparisons with other state of art detectors in different scales.

Dialogue Pre-processing
We manually labeled over 12, 000 dialogues for training and evaluation. They are divided into 4 categories: positive, nega-tive, wondering, and neutral. Most of them come from scripts of TV shows, movies and books. Numbers of dialogues in different categories are listed in Tab.1

Typical Routine
Before diving into the training process, normally several pre-processing steps should be done in order to improve overall efficiency. For English speech classification work, a typical routine is as follows: • Spell check: A widely used library for spell check is PyEnchant. This step is optional.
• Stemming or Lemmatization [6]: Both of them are used to find the original form of words, e.g. 'talk' is the original form of 'talking' or 'talked'.
• Case normalization: Convert input to uniform uppercase or lowercase. For example, "Today I'm happy." might be converted to "today i'm happy.". However, for Chinese classification or recognition system, a special step, segmentation, generally needs to be done. Normally researchers will build a LUT (look up table) which contains thousands of common used Chinese words before training, and then the frequency of appearance of each word in the target dialogue will be checked. A number of segmentation algorithms have been developed over the last decade [7][8][9]. The frequency vector is finally used as training materials.
One of the most appealing advantage with this method is that outputs of segmentation are of uniform length (equals to the length of pre-built LUT), so that during most classification algorithms, e.g. KNN and SVM, it can achieve a better results. Another important advantage is that segmentation with pre-selected keywords saves much time for reseachers and im-proves the efficiency [10].
However, drawbacks with this algorithm shouldn't be ignored.
• Over segmentation: Similar to image precessing, over segmentation may significantly decrease the accuracy of our classifier [11]. For example, 'not happy' may be segmented to 'not' and 'happy'. Obviously a single 'not' is meaningless but a 'happy' leads to a positive emotion, thus this sentence will be classified as positive because of over-segmentation, while the original emotion should be negative.
• Loss of linguistic features: Order of words is lost after segmentation. Same set of words with different orders sometimes gives absolutely opposite emotions. For example, as demonstrated in Fig.1, "He's happy while I'm sad.", and "I'm happy while he's sad" share a same vocabulary, but give opposite emotion.

Pre-pocessing in EmoNet
In order to avoid the problems listed above, we didn't adopt segmentation in EmoNet. And because our net is tested in Chinese dialogues, stemming and lemmatization are not necessary.

Removal of Stop Words
Many words or expressions do not contribute to emotion or have same like-hood of occuring in those sentences not relevant to the target one. Such words are called stop words [12]. In order to save storage space and improve classification efficiency [13], we need to remove them before further training. A list of Chinese stop words from Data Hall is used [14]. As shown in Fig.2, we searched for stop words in every source dialogue and removed them in this step. Notice that the order is kept.

Re-encode Input Characters
UTF-8 (8-bit Unicode Transformation Format) is a large and powerful set of characters which encompasses most of the world's writing systems [15]. With UTF-8, EmoNet gains the ability to perform in multi-language dialogue systems. Normally, one UTF-8 character is encoded with 1-6 byte, and in Chinese encoding, 3-byte UTF-8 is used.
However, because of UTF-8's large range, most information in this encoding system is useless. And because of the great number of Chinese characters, it's unnecessary to express each character precisely. The next step is to re-sample the encoding of input characters.
Firstly, because our net currently only contains Chinese, English and Arabic numerals, we can truncate UTF-8 to a new set that only contains characters of these vocabularies. Tab.2 gives encoding of Chinese characters, full-width English characters and full-width Arabic numerals in UTF-8. After this procedure, there are totally 20964 characters left, most of which are Chinese characters.

Encoding and Re-sampling
Unlike image processing system, encoding of characters is not "continuous", which means adjacent codes not necessary give similar meanings, while value in imaging systems gives continuous effect in hue, saturation, etc. Thus we don't need to represent each character precisely. In EmoNet, we re-sample the 20964 codes into 256 codes with Eq. 1.
After this procedure, every character is represented with one Byte. We set the maximum length of one dialogue to be 144 in EmoNet, thus totally 144 × 1 = 144 Bytes are used to represent one dialogue.

Convolutional Neural Network in EmoNet
We adopted a deep CNN to do the classification job in EmoNet. CNN is a powerful tool with an adjustable capacity by adjusting its depth and breadth. With large datasets, CNN can achieve high performance in learning the nature of speech systems [16].
In EmoNet, we adopted a CNN which shares a similar but deeper configuration with LeNet-5 [17], which consists of steps of convolutions, pooling, and fully-connected layers. Fig.3 gives the overall structure of EmoNet. Input layer is a 144 × 1 vector. The augmentation layer is used for feature expansion, which outputs a matrix of 32 × 32 × 3. Five convolutional layers are used, and most of them come with a max pooling layer for sub-sampling and dimension reduction. The last convolutional layer gives an output of 6 × 6 × 256. A fully-connected layer with dropout and softmax method is used for classification. The Output is of 4-dimension, positive, negative, wondering, and neutral.
Following sections give introduction of function and configuration of these layers in detail.

Augmentation Layer
We applied an augmentation layer to expand feature from the input and scale for the following layers. Because generally every convolutional layer outputs a matrix of smaller size but larger dimension compared with its input matrix, in order to apply enough number of convolutional layers, we need to expand the input matrix to a larger square matrix.
In every convolutional layer, we chose to use a filter of size 5 × 5 and stride 1. In order to outputs a 12 × 12 × n in the last step of convolutional step and apply 5 convolutional layers, we need an square output of size 32 in the augmentation step, which can be calculated with Eq. 2.

Convolutional and Pooling Layer
The convolutional part consists of 5 convolutional layers and 4 pooling layers.

Convolutional Layer
The convolutional layer adopts a filter to compute convolutions and generates a higher dimension matrix. In EmoNet, the filter sizes of all the 5 convolutional layers are chosen as 5 × 5. Filter dimensions vary in different layers.
Typically, a sigmoid nonlinearity (or activation function), which is given in Eq.4, is used after every convolutional layer.

Pooling Layer
Pooling layers introduce invariance, reduce dimension and prevent overfitting in CNNs [18]. We adopted overlapping pooling layers in this network. Overlapping pooling means adjacent pooling regions overlap with each other. For example, the pooling region size is denoted as s × s, and stride is denoted as z. If s ≥ z, which is used in most CNNs, it is non-overlapping. If s < z, then it becomes an overlapping pooling. According to Krizhevsky, Sutskever, and Hinton, an overlapping scheme will reduce error rate and decrease the probability of overfitting [20]. We adopted s = 5, z = 1 in EmoNet.
Also, many pooling methods have been developed, such as max pooling, average pooling, chunk pooling, etc. We used max pooling in EmoNet. Its principle is illustrated in Fig.4.

Fully Connected Neural Network with Dropout 3.3.1 Fully Connected Layers
A 3-layer fully connected neural network is connected to the last convolutional layer. Fully connected layers are affine layers, and every neuron in one layer is connected to each neuron in its adjacent layer. Each neuron of the first two layers comes with an RELU activation. Detailed configuration is listed in Tab.3. The final layer outputs a 5x1 vector. This vector is then transformed into 5 categories, positive, negative, wondering, neutral and meaningless (which means this input should be discarded) by softmax classification.

Overfitting
Overfitting is a common problem in neural networks, especially when researchers don't have a large enough dataset. Nowadays there are two popular methods to prevent overfitting, batch normalization and dropout.
We tested EmoNet with batch normalization first. During training, inputs of every layer are always changing with the fact that parameters in the prior layers are changing at the same time. This phenomenon, which is also called internal covariate shift [21], slows down the training process and requires careful initialization. Batch normalization solves this problem by normalizing layer inputs. However, in EmoNet, it reduced the error rate only by 2.5% to 3.8% compared with EmoNet without batch normalization with same amount of training. We think because our input is characters, which are discrete (while input of imaging systems are continuous), normalizing them might be unhelpful. Thus we tried dropout instead.
The core principle of dropout is to randomly disable some neurons, along with its connections during trainning [22], but do nothing during test. According to Srivastava et al., the optimal retention probability p of input layer should be close to 100% while with hidden layers it should be close to 50% [22]. We set p input as 1.0 and p hidden as 0.3, and we saw an reduction of error rate of around 6.5%.

Training Process
Parameters need to be updated during training, and normally people use back-propagation [24]. During this procedure, many parameters updating methods have been developed, such as SGD, momentum, RMSprop, and Adam. Adam is a first-order gradient-based optimization based on adaptive estimates of lower-order moments which can achieve a faster convergence than other methods like SGD, momentum and RMSprop [25].
We chose Adam as EmoNet's parameter updating method. We divided the dataset into 2 parts, one of which was used for training, the other was used for evaluation. Cross validation [26] is used during the parameter updating process, and the training data was further divided into 32 mini-batch for cross validation.

Test of Different Configurations
Configuration of EmoNet is modified from VGGNet [23]. We tested different configurations of EmoNet. Validation error and training speed within one epoch using same amount of data are used as evaluation criteria. Tab.4 gives the configurations we tested, and Fig.5 gives comparison between them. Matplotlib [27], a python plotting package, is used for drawing. We chose configuration B, which achieved relatively high accuracy and fast performance at the same time, as the final structure of EmoNet.

Test of Different Parameters
Learning rate, regulation strength in affine layers, and mean of initialization are also important hyper-parameters that give huge influence on final results. Because dropout makes careful initialization less important [22], we only tested different sets of learning rate γ and regulation strength L. Results are collected with same amount of input data and same number of training iterations. Results are shown in Fig.6.
The group of highest accuracy, group C is chosen, where γ = 5 × 10 −6 , L = 1.5 × 10 −4 . The loss v.s. time curve of the first epoch is shown in Fig.7.

Evaluation
With optimized parameters, EmoNet was trained with tensorflow and Nvidia GPU. After 100 epochs' training, an overall top-1 accuracy of 72.8% was achieved. Tab.5 gives accuracy of each category separately.
According to the results, neutral and wondering states are relatively easier to detect than neutral and negative states.
Wondering states are the easiest to detect, which is also easy to interpret. Generally dialogues of this emotion state come with question marks, or some specific words, like "what", "why", "how", etc.  Although accuracy of neutral state is high, we found that most emotion states of false detected dialogues were detected as neutral. This phenomenon may result from two facts, the first one is that a neutral state is too vague and general, and the other one is that unlike other states in which we can find specific emotional words (for example, "why" for wondering states and "happy" for positive states), there are no specific words for neutral states.
Detection with positive states and negative states are similar. However, as we can see, accuracy with positive states is low. This may result from the quality of the training materials. Ironies are very common in our source, and it is sometimes even hard for us to distinguish whether a dialogue is positive or negative.
Comparisons between EmoNet and state of art Chinese textbased emotion detectors are given in Tab.6, Tab.7, and Tab.8.
Multi-Model Net was given by Ze-Jing Chuang and Chung-Hsien Wu [28], and ESiN was given by Jianhua Tao [29].

Conclusion
In this paper, EmoNet was presented, an emotion detection system based on deep convolutional neural networks. We analyzed the feature of Chinese encoding and adopted a preprocessing step without segmentation, stemming or lemmatization, which introduces difficulties but address the problem of loss of linguistic features. Now in EmoNet, a simple resampling step was used to replace these steps. In the future, we will try some other algorithms. We are now working in developing a linear mapping system which can map arbitrarylength dialogues into equal-length outputs.
Different CNN configurations and different hyperparameters were tested. An overall accuracy of 0.72 was achieved with 12,000 training dialogues, 100 epochs' training and the optimized modal. Top-1 accuracy of EmoNet is higher than other Chinese text based emotion detector. However, it's foreseeable that with more training materials, EmoNet has the capacity to achieve a better performance.


extfieldsofstudy:
Computer Science
Medicine

text:
Predicting the bounds of large chaotic systems using low-dimensional manifolds

Predicting extrema of chaotic systems in high-dimensional phase space remains a challenge. Methods, which give extrema that are valid in the long term, have thus far been restricted to models of only a few variables. Here, a method is presented which treats extrema of chaotic systems as belonging to discretised manifolds of low dimension (low-D) embedded in high-dimensional (high-D) phase space. As a central feature, the method exploits that strange attractor dimension is generally much smaller than parent system phase space dimension. This is important, since the computational cost associated with discretised manifolds depends exponentially on their dimension. Thus, systems that would otherwise be associated with tremendous computational challenges, can be tackled on a laptop. As a test, bounding manifolds are calculated for high-D modifications of the canonical Duffing system. Parameters can be set such that the bounding manifold displays harmonic behaviour even if the underlying system is chaotic. Thus, solving for one post-transient forcing cycle of the bounding manifold predicts the extrema of the underlying chaotic problem indefinitely.

Introduction
Coupled systems of non-linear ordinary differential equations (ODEs) are ubiquitous in scientific literature and appear as the end result of efforts to mathematically model phenomena from practically every branch of science. The desire for quantitative accuracy often leads to such models containing large numbers of variables, for example, as a consequence of a highly resolved finite element mesh, an in-silico model of human tissue-be it heart, tumour or brain -containing many simulated cells, or an electrical circuit with many nodes. Systems of ordinary differential equations with non-linear terms have the potential to display chaotic behaviour, in which case their solutions evolve on "strange" attractors of fractal dimension.
The ability to precisely predict the extrema of large chaotic systems has obvious utility across many fields. While bifurcation analysis serves to characterise the qualitative behaviour of a system, extrema prediction adds quantitative knowledge. Perhaps the clearest benefit is where state variables have critical values, e.g., in predator-prey systems where reaching a a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 possible future states of a system relatively narrowly in phase-space and that they can be applied to high-D chaos in high-D phase space. The obvious drawback is their local and consequently short-term nature. The effort towards determining quantitatively accurate long-term global bounds of chaotic systems has mainly been focused on obtaining closed-form solutions to specific systems on an ad-hoc basis. As the information required to store a manifold is exponentially dependent on dimension, existing methods of bounding strange attractors break down in high-D phase space. Efforts to bound chaotic systems have, so far, mainly focused on systems of three variables [40,41] with attractor dimension close to phase space dimension. A more generalised method for the bounding of low-D chaotic systems was proposed in [42,43], which took the approach of defining localising sets with bounding manifolds defined as polynomials.
The present method is applicable in situations when the phase-space is high-D, but the attractor is low-D. In this situation it is computationally feasible to define a long-term bound for the system. Although study of high-D chaos and the factors impacting on attractor dimension has taken place since the 1980s [39,[44][45][46], as discussed above, many natural processes, even if high-D in terms of phase space, generate low-D attractors. It has been and still is commonplace to reduce large systems to the minimum dimension necessary for qualitative accuracy [5], which in turn tends to be defined as preservation of the topology of the attractor [9]. However, qualitative methods sacrifice quantitative accuracy. And quantitative accuracy may be crucial. Even periodic systems, are often simulated by very large quantitatively accurate high-D models. The study of dynamic phenomena with linear finite element models is an example of this.
The method presented here finds its niche in situations where model reduction is unacceptable and the extrema of system variables are important. It is deliberately designed to be intuitive and easy to apply, i.e. "wrap around" an existing system. The aim is that, once applied and proven on attractors of higher (but still low) dimension, it can be used by non-experts in chaos and non-linear dynamics in a plug and play fashion.

Results
Two manifolds are defined in the following: The auxiliary manifold (AM) and the bounding manifold (BM). The AM is an M-dimensional manifold, where M is an integer. It is embedded in the phase space of the system to be bounded and, in turn, the BM is embedded on the AM. The hope is that the AM and the BM display stationary harmonic behaviour in time, such that they only need to be calculated for a couple of post-transient forcing cycles.

The bounding-and auxiliary manifolds
The manifolds are non-autonomous extensions of the concept of invariant and inertial manifolds and draw inspiration from existing literature [47,48]. However, the following does not adhere strictly to the quite specific definitions and nomenclature of inertial manifold theory. Firstly, inertial manifolds are most often employed as a means to directly solve the original equations on the manifold. They allow the reduction of dimensionality from either infinity (for partial differential equations) or high-finite dimension (ODEs with many degrees of freedom). Here, no attempt is made at solving the underlying system of ODEs on the manifolds. The AM serves the purpose of reducing the dimensionality of the BM-not that of the underlying problem. Secondly, inertial manifold theory rests on the premise that solutions in full phase-space approach a solution on the manifold. The requirement here is simply that solutions approach the manifold, which is less restrictive. So, although the AM draws from inertial manifold theory, it is named as it is to emphasize the differences and avoid confusion.
If all unstable directions are contained on the AM, solutions that do not start out on the manifold, will decay exponentially onto it and, post-transiently, an attractor will lie on the AM. Let c denote the bounding manifold. Then, the premise of the method is that where A is an attractor, B is the auxiliary manifold and C is an M-manifold with boundary. Thus, c is a closed (M − 1)-manifold. Note that "bounding manifold" is formally a (potential) misnomer as extrema of the system may lie on interior points (and not on the BM itself). Further, this necessitates the post hoc calculation of interior points, which is discussed below.
The AM is embedded in the full phase space x 2 R N of the system of ODEs: As it shall become apparent in the following, it is a requirement that f is smooth. Furthermore, it is a requirement that M > d, where d is the (fractal) dimension of the attractor. The most obvious choice is to select M as the nearest integer greater than d. As mentioned above, the discipline of bounding chaotic systems has previously focused on N − 1 dimensional boundaries to systems with N dimensional phase spaces [40,42,43]. Such approaches will fail rapidly as N increases to that of even modest sized problems, since the computational cost of operating on a discretised manifold depends exponentially on its dimension. Here, it is exploited that the strange attractor is normally a low-D entity in a high-D phase space, i.e, that d < M < <N. So, the BM, c(t), need not be more than M − 1 dimensional. If c(t) is harmonic in time, the solution can be stopped after a couple of post-transient forcing cycles. Fig 1 shows this scenario for M = 2. For the method to work, the AM must be invariant under the flow, meaning that solutions that start on it, stay there. The manifold is allowed to vary with time. So, labelling the manifold "invariant" may cause confusion. Nevertheless, this is the naming convention widely adopted in literature and time dependent inertial manifolds are well known [49]. The key requirement is that the manifold displays regular motion. Clearly, if the manifold itself is chaotic, hence the BM also, the method has achieved very little.
To introduce key concepts, we consider the case of a 2-D AM in a R N phase space and then generalise the ideas to M-D manifolds in R N phase spaces. Fig 2 shows a schematic of the BM and AM. It also introduces the basis vectors of the AM (w 1 ,w 2 ), the inward normal u 1 and the BM unit tangent vector u 2 . Furthermore, the figure shows an off-manifold point which is projected along a stable direction onto the manifold. If certain criteria are met, it can be proven that the off-manifold point converges exponentially with a solution on the manifold [9,48,50,51]. However, what matters is that off-manifold points approach the manifold asymptotically, end up somewhere on it and stay there, which is uncontroversial. Here, the notion of points being "on" the manifold is taken to include points infinitesimally close to it, which is the case with asymptotically approaching points. Obviously, once the distance drops below machine precision, there is no practical distinction.
The motion of the BM is the sum of the tangent on-manifold motion and the orthogonal off-manifold motion. The latter is readily determined through Gram-Schmidt processing of the flow. The expression for BM motion then becomes: where f(x, t) = f(c(m), t) is defined in Eq (2), _ cðm; tÞ is the rate of the BM, d i are the on-manifold velocities and u i is a local orthonormal basis spanning the tangent space of the AM. For AMs of dimension 3, 4 and so on, Eq (3) is easily expanded with d 3 , u 3 , d 4 , u 4 , etc. That is, generally The flow and the vectors are evaluated locally along the length of the BM and depend on time, i.e. Initially they lie as a cloud in phase space, but they rapidly contract along stable directions until they converge to the AM. As shall become apparent, it is sufficient to calculate the tangent space of the AM to orientate the BM c(m, t). This is done with a spectral method.
The fundamental idea behind the method is that the locally unstable directions, the directions of stretch, are contained in the on-manifold terms ("k"). If that is the case, the remaining offmanifold motion ("?") will consist only of exponential decay onto the manifold.
The nomenclature is such that u 1 is the inward normal (tangent to the AM and orthogonal to the BM) with all other vectors spanning the space tangent to the BM (tangent to both the AM and the BM). The inward motion of the boundary is taken to follow the heuristic expression: where β is the linear impedance and κ = κ(c(m, t)) is the signed dimensionless curvature signal passed through a low-pass filter with time constant τ: The equation of motion Eq (6) automatically bounds the underlying ODE and filters away high frequency oscillations. It will straighten the BM locally until the curvature term vanishes, at which point the BM is at equilibrium. Since the BM is a closed curve, the global net effect of the curvature term is to contract the BM. The responsiveness of the BM is determined by the values of β and τ. Although the expression is non-linear in d 1 , it has the merit of being local, so does not require the (iterative) solution of a system of non-linear equations. A spatially coupled formulation, e.g., based on constrained optimization and a Lagrange multiplier field, would come with this drawback. The present study considers 1-D BMs, so the measure of curvature is given. For higher-dimensional BMs (2-D surfaces and beyond), a scalar quantity such as the mean curvature might be used. The motion of the BM is discretised using Hermite finite elements. These have the sufficient order to allow the calculation of curvature. Axial motion. The axial motion of the BM seeks to space nodes evenly along the BM. This is no more than an adaptive meshing exercise. For the sufficiently well discretised BM, the axial motion has no appreciable impact on the shape of c. The below expressions are valid for a 1-D BM, so a general (M − 1)-D version is needed. Consider the rate of change of elastic energy with the addition of a dissipative term Now, applying Eq (3) and requiring stationarity with respect to the d 2 field, one gets the expression where prime denotes differentation with respect to m. As with the d 1 field described above, the d 2 field is expanded in terms of Hermite shape functions and the finite element method is used to solve Eq (9) for nodal values of d 2 . The details are given in the methods section.

Interior points
Unlike in some other methods from the literature, the interior of the manifold is not necessary for the determination of manifold orientation. Despite this, it must still be considered, as an extremum of a state variable might lie on the interior and not on the BM. One must balance the added computational cost of sampling the interior against the need for a high number of samples. The equation of motion for the interior points is a compromise between simplicity and the desire for evenly spaced points and could be defined in a variety of ways. What is unequivocal though, is that interior points must ride on the tangent space associated with the AM. The orientation at interior points is determined in the same way as on the BM. A practical and fairly straight forward way of ensuring evenly spaced points throughout a simulation, where the BM may change shape significantly, is to connect the interior points to their neighbours and the BM with linear truss elements of zero initial length. Then, point positions are updated dynamically as nodes of truss elements. For the purposes outlined here, this is nothing more than a simple, and probably sub-optimal, adaptive meshing technique. The advantage of truss elements is that they are readily defined in higher dimensions as their length can be easily calculated in N-dimensional Cartesian coordinates. Thus, they are practical for point spacing in phase spaces above 3-D. The equation of motion for interior points is identical in form to Eq (3), except that now d 1 and d 2 are determined by force summations of the truss elements. That is, the interior points, i.e., truss element nodes, are restricted to move on the AM onto which all truss element forces, accelerations and velocities are projected. The truss element normalised stiffness and damping, k t /m t and d t /m t , respectively, are given in Table 1. The motion of truss element nodes on the BM c is prescribed to follow the BM in a one-way coupling fashion, so the motion of the BM is not influenced by the interior nodes. This means that the motion of the BM can be computed first, and interior point sampling performed in postprocessing.
The equation of motion for an interior point, where d i are determined as and h 2 R N is the resultant mass-normalised force from trusses pulling on an interior point.

Example: The duffing oscillator
The Duffing oscillator is among the simplest systems known to exhibit chaos, yet still finds use for a wide variety of applications including, e.g., nano-mechanical systems [52]. Consider a system consisting of a Duffing oscillator in series with a number of linear spring-mass-damper oscillators. The linear oscillators add to the phase space of the system, and quantitatively affect Table 1. Model parameters for the simulations. The parameter τ is the time constant for curvature filtering, β is the proportional impedance on curvature, β 2 is the impedance on axial motion, k t , d t and m t are the stiffness, damping and mass of interior point truss elements, α is the non-linear Duffing stiffness, k 1 , d 1 and m 1 are the linear stiffness, damping and mass of the Duffing oscillator, f 0 is the forcing term amplitude, ω is the forcing term frequency, k 2 , d 2 and m 2 are the stiffness, damping and mass of the linear oscillators, Δt is the time step for the time integration scheme (BM and interior point position update), Δt λ is the time step for eigenvalue calculation and orientation update, K is the number of nodes (and elements) used to discretise c, and K G is the number of Gauss points for numerical integration in space along c.
All the solution to the equations, but do not add to the dimension of the attractor or influence its topology. If the duffing oscillator is assigned the first two degrees of freedom, the system of equations can be written as where f 0 cos(ωt) is a harmonic forcing term, k 1 , d 1 , α and m 1 are the parameters of the Duffing oscillator and all the linear oscillators are assigned the identical parameters k 2 , d 2 and m 2 . The governing equations are on the form of Eq (2) so can be plugged into the method as described above. In order to reduce the mesh refinement and time step requirements associated with curvature and stretch along fast, linear terms, the equations are transformed as ., x N ] T = [y 1 , y 2 , y 3 /10, . . ., y N /10] T and transformed back in post-processing. This transformation is beneficial and possible because it is known a-priory that the non-linearity of the system appears in the first two equations-thus "compressing" all other variables and effectively increasing mesh refinement in the non-linear directions. For some systems, all variables might appear in non-linear terms, or it may not be known which ones that do. Then, scaling becomes futile, potentially necessitating finer meshes and shorter time steps. Parameters for all simulations are shown in Table 1. In general terms, a fast system requires a fast BM. So a system containing high frequencies and rapid changes in J will require low impedances on the BM, lower time constants for the first order filter and, of course, smaller time steps. In Fig 3 is shown the convergence of the AM with twenty solutions for parameter set II in Table 1. The phase space is 42-D, and only the (y 1 , y 2 , y 42 ) projection is shown, so more information is contained in the simulation than can be displayed on the figure. Note that the solutions converge to the AM and stay there. The figure also shows the truss elements connecting interior points to each other and to the BM. A video of a similar process for parameter set III is shown in S1 File. Fig 4 shows snapshots of projections of the BM onto a selection of coordinate surfaces for parameter set III in Table 1. The projections vary considerably in shape and size, indicating differences in phase and magnitude, respectively, of the harmonic components of the state variable motion. A video of a similar process for parameter set III is shown in S2 File.  Fig 3). The motion of the AM and the BM is calculated for five forcing cycles and the subsequent 95 cycles are copies of the fifth, which is post-transient. Twenty time series are calculated for 100 forcing cycles in order to validate the BM and it is seen that they all remain within the predicted values. The BM is noted to be reasonably tight and to be following the low frequency dynamics of the system. As one would expect, not all variables in the system influence the BM equally. The geometric nature of the method means that for the BM to contract, curvature must exist in the coordinate direction of the relevant variable. With the proposed equation of motion, the BM will tend towards being convex everywhere. The approach works well with the here presented Duffing-type system, but there is no proof that it will in general, so future modifications to the equation of motion Eq (6) producing a locally concave BM might be beneficial.

Example: The parametric duffing oscillator
To illustrate the method's ability to deal with situations where the Jacobi matrix depends explicitly on time, a parametrically excited version of the duffing equation is considered. The equations for this system are given in Eq (13). Notice that the non-linear term now varies explicitly with time.

Local and instantaneous stretch and compression
Let us explore the behaviour of two infinitesimally nearby trajectories in a smooth, non-linear and non-autonomous flow. Due to smoothness, the state can be represented by its Taylor expansion. Consider a perturbation in time by Δt and in space by a v, where a is a scalar and v 2 R N is a unit vector: where Δt = t − t 0 and O(3) represents terms of order 3 and above. Fig 10 shows a state, x 1 , perturbed forward in time and another, x 2 , perturbed in space and time. Representing the perturbed states x 1 and x 2 as their Taylor expansions and substituting _ x ¼ f we obtain where f = f(x 0 , t 0 ). Now, from Eq (15) and recognising that @x @a ¼ v, the rate of stretch is readily obtained as lim a!0 It is easily seen that Where J is the Jacobi matrix. If v is an eigenvector with a real eigenvalue, the local stretch is parallel to it. If it is complex, the stretch will lie in the space spanned by v and its complex conjugate (since the flow is real). Thus, imagnining for a moment a full spectral decomposition, it is easily seen that, locally, the only directions of positive stretch are those associated with eigenvalues with a positive real part. Nearby solutions will converge exponentially in all other directions-i.e. contract. This holds for either autonomous or non-autonomous systems as no assumptions have been made in regards to the time dependence of _ xðx; tÞ apart from the requirement of smoothness. So, as is evident, the spectrum plays a pivotal role in manifold orientation, even for non-autonomous systems.
Although the above is not a novel result, its implications justify its repetition here. Chaos is a global phenomenon and eigenvalue decompositions cannot be utilised globally to describe stretching and folding. However, locally and instantaneously, stretching and folding are described by the spectrum of the (ever changing) Jacobi matrix and the recognition of this underpins the present method. Spatial discretisation. The choice of spatial discretisation method is not central to the method, and the approach taken here may well be substituted with another method depending on application. Here, all fields on the BM are discretised using finite element theory with Hermite shape functions: where ϕ i (m) are the shape functions, K is the number of degrees of freedom used to discretise the BM and, e.g., are nodal values of the position of the BM and N is the dimension of the system. Time derivatives and spatial derivatives for the BM are easily obtained from these expressions. The use of Hermite (higher order) shape functions allows for the definition of curvature, which appears in Eq (6). All simulations use the same 40-node (K = 80) mesh of 1-D elements with isoparametric Hermite shape functions. Numerical integration over m is performed with a full Gauss integration scheme with 7 Gauss points per element (K G = 280). Nodal updates are determined by solving for nodal rates and integrating numerically over t with an explicit time integration scheme. The forward Euler method is used for ease of implementation. Introducing the approximations in Eq (18), multiplying Eq (3) by ϕ j and integrating over m on both sides, we obtain I Which can be written in matrix form and solved as where Now, the time integration scheme can be utilised to update c. The system matrices from the finite element formulation do not change between time steps, which can be exploited for efficiency.
For a higher order attractor, the expression for the right hand side is modified according to the attractor dimension M as The inward velocity d 1 can be evaluated at each point on c directly using Eq (6). The tangential velocity d 2 is found by the finite element solution to Eq (9): where Auxiliary manifold orientation -concepts. There exists an extensive body of literature on how to define and practically determine invariant manifolds. A comprehensive overview is given by Gorbana and co-workers in [53]. This is an active area of research and new highly advanced methods are continually being developed [54]. Here, an approach is taken, which is somewhat related to what is labelled the "method of invariant grids" in [53]. It differs substantially in that the tangent field of the invariant manifold is not found by numerical differentiation in phase space. This approach, though intuitive, requires refined meshes and is prone to numerical instability. Instead, a hybrid approach is taken, where AM orientation is determined by a spectral method, but selection of eigenvectors from this spectrum, and thus orientation, makes use of the grid. The reasoning behind this choice is pragmatic and based on experience: Even if very large time steps are taken between two orientation updates, the tangential motion of the AM in between these updates will be small when compared to typical mesh refinement. Furthermore, the use of a spectral method renders the orientation of the inward normal u 1 independent of interior points. The obvious drawback is the need for solution of the eigenvalue problem at each Gauss point and at every time step. A redeeming feature is that this can be done efficiently by use of sparse iterative solvers that solve for only the M necessary eigenvalues and utilise eigenvalues from previous time steps as initial guesses. The current software implementation does not exploit these features and computes the entire spectrum from full matrices using direct methods.
Locally, the tangent space of the AM is spanned by the tangent vectors w i 2 R N where i = 1, 2, . . ., M. Chaotic motion is often described as successive processes of stretching and folding. As discussed, locally, these processes correspond to exponential growth and decay, respectively. The selection of the tangent vectors w i is such that every direction in which stretch occurs is contained within the AM. In this way, every off-manifold motion will be an exponential decay onto the manifold.
Consider the Jacobi matrix of Eq (2) In any non-linear system, the Jacobi matrix depends on the state, i.e., J = J(x). In non-linear non-autonomous systems it may also depend explicitly on time, i.e. J = J(x, t). The method is equally applicable in both of these situations. The AM is oriented such that all directions of stretch are contained within it at all times. Locally and instantaneously, these directions are those spanned by eigenvectors associated with eigenvalues with a positive real part. That is, the manifold tangent vectors w i are selected such that they span the same space as solutions to the local system associated with directions of stretch. The number of eigenvalues with positive real parts determines the required dimension of the AM. For two positive real parts at least a surface is required, and for three a 3-D hypersurface, and so on. Knowing the dimension of the attractor a priori is advisable, since then a good (if not perfect) indication of the dimension of the AM is then known. A good account of the estimation of attractor dimension, for both experimental data and systems with known ODEs, is given by Wolf and co-workers [55].
When a direction of stretch with associated eigenvalue λ is encountered somewhere in phase space, there are two situations where v 1 and v 2 are the eigenvectors associated with eigenvalues λ and " l.
2. l 2 R: The direction of stretch is parallel with the eigenvector associated with λ.
In case 1, one may define the local tangent vectors in any way that span the same space as v 1 e lt þ v 2 e " lt , e.g., by selecting two distinct phases: 0 (giving v 1 + v 2 ) and +/−i (giving i v 1 − i v 2 ). If the flow is real v 1 and v 2 appear as complex conjugate, so the space is spanned by Re(v 1 ) and Im(v 1 ). These vectors are then Gram-Schmidt orthonormalised to yield the orthonormal basis (w 1 ,w 2 ) of the AM. With the requirement of smoothness, the orientation of the AM is defined even in regions of phase space with no directions of stretch. This implies the requirement that J(x) is continuous and differentiable. Even if a direction of stretch is associated with a positive real eigenvalue, there may be another region of phase space where this eigenvalue becomes complex (with positive or negative real part) and thus requires a second direction of expansion. The only way to know this is to trace the eigenvalue along the AM and track its evolution. A schematic illustrating this process is shown in Fig 11. The phenomenon of directions of stretch varying in number across phase space is well known and is what is labelled "unstable dimension variability" in, e.g., [39].
Auxiliary manifold orientation-initialisation. Selecting the appropriate eigenvectors during initialisation takes a two step approach: First, the eigenvalues of all Gauss points on the BM are determined until a node is reached where at least one eigenvalue has a positive real part. Subsequently, the rest of the Gauss points are scanned incrementally along the BM selecting the ones corresponding with minimum increments in modulus. If the discretisation of the BM and AM is sufficiently fine, the minimum increment selection criterion will correspond to continuity. Thus, eigenvalues for neighbouring points on the AM and the BM can be directly organised and numbered consistently. Although it is somewhat computationally wasteful, it is conceptually more straight forward to perform the sorting process before the selection of eigenvalues. This is illustrated as pseudo-code in algorithm 1. Now, the number of entries in l j equalling 1, should exceed the a-priori predicted dimension of the attractor. If it does not, then the initialisation of the BM is inadequate and another must be tried.
The above pseudocode is written to communicate the concept of the method and is not optimised. For example, there is no need to sort eigenvalues that do not correspond to directions of stretch.
Auxiliary manifold orientation-updates. Updating l ðpÞ j at each time step simply follows the principle of continuity in time, so updates in eigenvalues are those minimizing changes in modulus from the previous time step. At each time step, the inward normal u 1 is projected onto the orthonormal basis associated with the corresponding eigenvectors. Thus, at each time step, Gauss point and interior point the Jacobi matrix is updated, the eigenvalues of the Jacobi matrix are calculated, and the eigenvectors associated with minimum change in eigenvalue modulus are used to update the orientation of the manifold. If Δt λ is sufficiently small, this corresponds to continuity in time of manifold orientation. The process of time integration and manifold updates is outlined as pseudo-code in algorithm 3. Note that the time step between orientation updates is larger than the time step for position updates. For clarity, the algorithm is written for a 2-D AM (M = 2) with l ðpÞ 1 and l ðpÞ 2 corresponding to the eigenvalues of orientation at gauss point m p . Generalisation to higher dimension is straight forward. J = J(c(m), t i ) Solve J v j = λ j v j for λ j and v j , j = 1, 2, ..N Select λ j from the spectrum of eigenvalues such that jl ðpÞ 1 À l j j is minimised l ðpÞ 1 l j Select λ k from the spectrum of eigenvalues such that jl ðpÞ At the first time step, updates in orientation are associated with large changes in u 1 if the initial guess is poor. Conversely, if the initial guess for the AM is close to the real, long term, AM, the projection will not alter u 1 much. In this paper, a flat, circular surface is used as an initial guess and initial values of w 1 and w 2 computed on this surface. This approach works well when off-manifold directions are almost in phase with on-manifold location, but may fail if off-manifold dynamics are out of phase, potentially causing the AM to become degenerate and numerically unstable. Thus, convergence depends on a reasonable initial guess. This which might potentially be generated from Poincaré sections in future incarnations of the method. Fig 3. shows a manifold changing its shape in the initial transient phase of a simulation. Fig 12 depicts the converged manifold and the manifold eigenvalues for parameter set III in Table 1. The manifold is projected to the (y 1 , y 2 , y 81 ) space. Where complex, as expected, the two eigenvalues of orientation are complex conjugate. There is a central region of stretch where λ 1 is positive. Conceptually, this corresponds to the area around the unstable saddle equilibrium at (0, 0) in the unmodified 2-D Duffing system.

Discussion
A practical method for calculating extrema of chaotic systems has been presented here. The equations of the method are written such that they can be wrapped around an underlying system of which the detailed structure is less important, as long as its flow is continuous and differentiable. The simulation results are post-transient extrema of systems with fairly fast, close to in-phase, off-manifold dynamics. Such systems pose the least challenging in terms of defining (guessing) an initial shape, that will rapidly converge onto the real AM. Addressing this challenge, for example by use of fits to Poincaré sections, would increase the chances of convergence.
The presented method has drawbacks which should be addressed in future work: Firstly, it will not distinguish between basins of attraction as the BM is defined a priori as a single closed curve. For this, an approach based on the level set method would be superior but perhaps at the cost of added computational work. Secondly, the presence of noise would also make the BM itself noisy. The method in its current incarnation offers no strategies to deal with this added challenge, which will inevitably come with the application to experimental data. Thirdly, the adaptation to experimental data would entail a phase space reconstruction step. The method, in its current incarnation, requires knowledge of the rates of all phase space variables, so will need to be modified if phase space reconstruction techniques based on time-delay coordinates [56,57] are to be employed. However, it may well be that a hybrid scheme where some rates are inferred from knowledge of the rest would suffice.
The computational benefits associated with defining the BM as low-D as possible are potentially enormous. Consider the cost of operating on a conventional N − 1 BM rather than an M − 1 ditto. The memory needed to store the BM, let alone operate on it, scales as K D , where D  (12). The parameters (set III) are given in Table 1.
https://doi.org/10.1371/journal.pone.0179507.g012 is BM dimension and K is the resolution of the BM, e.g., as in this case, the number of nodes of a finite element formulation. Plugging in the numbers for even modest sized systems quickly reveals why the D = N − 1 approach is infeasible. For example, storing an N − 1 BM, for an N = 82 phase space system with identical resolution to the one used in the presented analyses requires 1.13 Á 10 155 bytes if stored as double precision.

Conclusion
Based on the presented results, it is concluded that it is practically possible to predict extrema for large differentiable chaotic systems and that these extrema have the potential to display harmonic behaviour. This enables the efficient calculation of extrema of chaotic systems in the long term once the post-transient regime is reached.