text:
Alexander N Pchelintsev
email: pchelintsev.an@yandex.ru

On a High-Precision Method for Studying Attractors of Dynamical Systems and Systems of Explosive Type

Keywords: attractor, high-precision calculations, power series, asymptote, population explosion model, piecewise smooth system, the Chua system, the number π, systems with hysteresis

The author of this article considers a numerical method that uses high-precision calculations to construct approximations to attractors of dynamical systems of chaotic type with a quadratic right-hand side, as well as to find the vertical asymptotes of solutions of systems of explosive type. A special case of such systems is the population explosion model. A theorem on the existence of asymptotes is proved. The extension of the numerical method for piecewise smooth systems is described using the Chua system as an example, as well as systems with hysteresis.

Introduction

Let us consider a dynamical system with quadratic nonlinearities of the form

Ẋ = B 0 + B 1 X + φ(X), (1)
where X(t) = [x 1 (t) . . . x n (t)] T is a vector function of time t with values in space R n , B 0 ∈ R n is a given column vector, φ(X) = [φ 1 (X) . . . φ n (X)] T , φ p (X) = ⟨Q p X, X⟩, B 1 and Q p (p = 1, . . . , n) are the matrices (n × n) of real numbers.

There has been a strong interest in the systems of the form [START_REF] Lorenz | Deterministic Nonperiodic Flow[END_REF] in scientific literature starting somewhere in the middle 20th century. The most popular object of study from the class of systems (1) is the Lorenz system [START_REF] Lorenz | Deterministic Nonperiodic Flow[END_REF]:

       ẋ1 = σ(x 2 -x 1 ), ẋ2 = rx 1 -x 2 -x 1 x 3 , ẋ3 = x 1 x 2 -bx 3 .
(2)

For this system, the matrices have the form:

B 0 = 0, B 1 =     -σ σ 0 r -1 0 0 0 -b     , Q 1 = 0, Q 2 =     0 0 -1 0 0 0 0 0 0     , Q 3 =     0 1 0 0 0 0 0 0 0     .
Such systems are of interest to researchers because they are systems with energy dissipation. There is compression of the volume of the phase space in them, and all trajectories are uniformly bounded. Thus, there is an attracting set called the attractor of the dynamical system to which all solutions tend at t → ∞. Consequently, the attractors of dissipative systems determine their behavior over long time intervals. However, for many systems of the form (1), the limit solutions are unstable. Therefore, correct numerical integration requires the high-precision methods that make it possible to flexibly control the accumulation of calculation errors over a long time interval. For example, with an increase of the order of the Runge-Kutta methods, the complexity of the formulas for calculating the approximate values of the phase coordinates (they have a tree structure) [START_REF] Hairer | Solving ordinary differential equations I: nonstiff problems[END_REF][START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF] grows.

Recently, to simulate of dynamic chaos in systems of the type [START_REF] Lorenz | Deterministic Nonperiodic Flow[END_REF], some researchers have begun to use effective modifications of the power series method. Gibbons in his article [START_REF] Gibbons | A Program for the Automatic Integration of Differential Equations Using the Method of Taylor Series[END_REF] proposed the quick formulas for calculating the coefficients of power series for the main types of right-hand sides of differential equations in the scalar form. Today this idea was generalized in a recursive procedure (called as automatic differentiation) to compute the values of the derivatives for power series [START_REF] Rall | Automatic differentiation: techniques and applications[END_REF]. An advantage over the general Taylor series method is that the calculations can be constructed by fast formulas in comparison to the direct symbolic differentiation of right-hand sides of nonlinear ODEs which requires a lot of computer memory for high-precision calculations. The method of power series in [START_REF] Hashim | Accuracy of the Adomian Decomposition Method Applied to the Lorenz System // Chaos[END_REF][START_REF] Abdulaziz | Further Accuracy Tests on Adomian Decomposition Method for Chaotic Systems // Chaos[END_REF][START_REF] Al-Sawalha | On Accuracy of Adomian Decomposition Method for Hyperchaotic Rössler System // Chaos[END_REF] is applied as the Adomian decomposition method (ADM). The Clean Numerical Simulation (CNS) [START_REF] Liao | On the Clean Numerical Simulation (CNS) of Chaotic Dynamic Systems[END_REF] is based on the method of power series at arbitrary-order and used the multiple-precision data, plus a check of solution by means of an additional computation using even smaller numerical noises.

There is the FGBFI numerical method (the firmly grounded backward-forward integration method) described in [START_REF] Lozi | A New Reliable Numerical Method for Computing Chaotic Solutions of Dynamical Systems: the Chen Attractor Case[END_REF][START_REF] Lozi | A New Accurate Numerical Method of Approximation of Chaotic Solutions of Dynamical Model Equations with Quadratic Nonlinearities // Chaos[END_REF][START_REF] Pchelintsev | An Accurate Numerical Method and Algorithm for Constructing Solutions of Chaotic Systems[END_REF][START_REF] Pchelintsev | On the Poisson Stability to Study a Fourth-Order Dynamical System with Quadratic Nonlinearities[END_REF]. In this articles, the recurrence relations for calculating of the coefficients of expansion of local solutions in a power series are received for any dynamic system with quadratic nonlinearities in the vector form (it distinguishes the FGBFI-method from automatic differentiation and CNS). This form of coefficients of the power series is simpler and faster to compute than in ADM (because it does not contain factorials). The author of this article derived the simple formulas of calculating the length of the integration step in the general form (it distinguishes the FGBFI-method from CNS). The criteria for checking the accuracy of the approximate chaotic solution are obtained. This article will consider this method which uses high-precision calculations to construct approximations to unstable solutions of the system (1), as well as its extension to piecewise smooth systems composed of quadratic nonlinearities (for example, the Chua's piecewise linear system).

It is also noteworthy that for certain types of matrices B 1 and Q p and the vector B 0 , the system (1) can have non-extendable solutions of the explosive type. For example, for scalar systems of the form

ẋ1 = 1 + x 2 1 (3)
with a particular solution x 1 (t) = tan t one of the asymptotes is

t as = π 2 . ( 4
)
The existing classical numerical methods available in the standard implementation of such mathematical packages including Mathcad, MATLAB, Maple, etc., cannot be used to build approximations to such solutions, since the numerical scheme will jump over the asymptote, "not feeling" it. But as will be shown in this article, the FGBFI method allows to get arbitrarily close to the asymptote without skipping the time t as , i.e., localizing it.

The FGBFI Method

The instability of the limiting solutions of the system (1), in fact, makes incorrect to use classical numerical methods over long time intervals. The limit sets of dynamical systems are constructed on such segments. We can use high-precision calculations but then we will have rather small integration steps. The solution to the problem is to apply a modification of the power series method with recurrent calculation of expansion coefficients.

Let us represent the solution of the system (1) as a series

X(t) = ∞ ∑ i=0 Υ i t i , (5)
where the vector Υ 0 is determined by the given initial condition. The remaining expansion coefficient vectors (5) are calculated by the following recursive formulas [START_REF] Lozi | A New Accurate Numerical Method of Approximation of Chaotic Solutions of Dynamical Model Equations with Quadratic Nonlinearities // Chaos[END_REF]:

Ψ i = [ψ i,1 . . . ψ i,n ] T , ψ i,p = i ∑ j=0 ⟨Q p Υ j , Υ i-j ⟩, p = 1, . . . , n, Υ 1 = B 0 + B 1 Υ 0 + Ψ 0 , Υ i = B 1 Υ i-1 + Ψ i-1 i for i = 2, 3, 4, . . . . (6)
The convergence region [-τ(Υ 0 ), τ(Υ 0 )] of the series ( 5) is limited for many nonlinear systems of the form (1), depends on the initial condition, and can be estimated by the following formulas [START_REF] Lozi | A New Accurate Numerical Method of Approximation of Chaotic Solutions of Dynamical Model Equations with Quadratic Nonlinearities // Chaos[END_REF]:

h 1 (Υ 0 ) = ∥Υ 0 ∥, µ = n • max p=1,...,n ∥Q p ∥, h 2 (Υ 0 ) = ∥B 0 ∥ + (∥B 1 ∥ + 2µ)h 1 (Υ 0 ) + µh 2 1 (Υ 0 ), if h 1 > 1, ∥B 0 ∥ + ∥B 1 ∥ + µ otherwise, τ(Υ 0 ) = 1 h 2 (Υ 0 ) + δ , t ∈ [-τ(Υ 0 ), τ(Υ 0 )], (7)
where δ is any positive number (can be chosen arbitrarily small). In these formulas, one of the compatible matrix norms ∥ • ∥ 1 or ∥ • ∥ ∞ with a vector norm is meant. The recursive formulas [START_REF] Hashim | Accuracy of the Adomian Decomposition Method Applied to the Lorenz System // Chaos[END_REF] for calculating the coefficients of the power series expansion were obtained from the substitution of power series into the system (1), as well from the products in Cauchy form. A theorem on the convergence of a power series is proved in the article [START_REF] Lozi | A New Accurate Numerical Method of Approximation of Chaotic Solutions of Dynamical Model Equations with Quadratic Nonlinearities // Chaos[END_REF]. The proof is based on the method of mathematical induction. The specified region from formulas ( 7) is obtained as a corollary from this proof.

Next, we describe an algorithm for constructing a trajectory arc on a given time interval [0, T].

In the papers [START_REF] Pchelintsev | An Accurate Numerical Method and Algorithm for Constructing Solutions of Chaotic Systems[END_REF][START_REF] Pchelintsev | On the Poisson Stability to Study a Fourth-Order Dynamical System with Quadratic Nonlinearities[END_REF], the numerical method based on the representations ( 5)-( 7) is called the FGBFI method. Next, we describe the algorithm underlying it.

Let T is the given length of the integration interval, S a is the ball containing the attractor of the system (1). Let us set such a representation of a real number that the machine epsilon

ε m << ε pw ,
where ε pw is an accuracy of evaluation of the common term of the series [START_REF] Rall | Automatic differentiation: techniques and applications[END_REF]. Thus, the summation by the Formula (5) stops at such a value i = i * when

∥Υ i * ∥ • |∆t| i * < ε pw , ( 8
)
where ∆t is an integration step. Note that for the convergence of the series, the value ∆t must be chosen as

0 < ∆t ≤ τ(Υ 0 ) in forward time or -τ(Υ 0 ) ≤ ∆t < 0 in backward time.
Since the convergence region of the series ( 5) is limited, we cannot take ∆t = T (as mentioned above, the value of T is large). Therefore, the trajectory arc on the time interval [0, T] will consist of the connected arcs for which the series (5) converges.

Since in the computer simulation we can only operate with a discrete set of phase coordinate vectors (instead of a continuous function X(t)), then, by analogy with the Runge-Kutta methods, we will obtain approximate values of phase coordinates through an integration step, only calculations will be carried out at each step using Formulas ( 5)-( 8) for t = ∆t = τ(Υ 0 ).

Note that here the error accumulation at each step is much smaller than for the fixed-order Runge-Kutta methods, since we can quickly enough calculate the missing terms in the sum using the Formulas ( 5) and ( 6) by reducing the value of ε pw . It can be seen from the formulas ( 7) that with such a numerical integration, the step will be variable.

Let the value way is the direction of integration: way = 1 for forward time, way = -1 for reverse time. Next, we show an algorithm for obtaining approximate values of phase coordinates.

Algorithm 1

The FGBFI Method 1: Set the vector X 0 ∈ S a of initial condition and the value way 2: ct := 0 ▷ The current time 3: ended := false ▷ Are we finishing the algorithm? 4: l := 1 ▷ Number of closed interval of series convergence (5) 5: Calculate the value of integration step ∆t = τ(X l-1 ) using the formulas [START_REF] Abdulaziz | Further Accuracy Tests on Adomian Decomposition Method for Chaotic Systems // Chaos[END_REF] 6: if ∆t > Tct then Calculate the vector Υ i using the formulas (6) 20: We moved beyond the ball S a

X l := X l + Υ i • p ▷ Add to X l

26:

Print "Decrease the value ε pw and/or ε m "

l := l + 1 36: GoTo Step 5
The criteria for checking the accuracy of the resulting approximate solution are described in [START_REF] Lozi | A New Reliable Numerical Method for Computing Chaotic Solutions of Dynamical Systems: the Chen Attractor Case[END_REF][START_REF] Lozi | A New Accurate Numerical Method of Approximation of Chaotic Solutions of Dynamical Model Equations with Quadratic Nonlinearities // Chaos[END_REF][START_REF] Pchelintsev | On the Poisson Stability to Study a Fourth-Order Dynamical System with Quadratic Nonlinearities[END_REF]. Note that the backward time pass is used in some of them.

The Software Implementation of the FGBFI Method and Numerical Simulations

For the software implementation of Algorithm 1 with

B 0 = 0,
and for the study of the limit points of dynamic systems of the form (1) for Poisson stability and the calculation of the Lyapunov exponents by the modified Benettin's algorithm [START_REF] Pchelintsev | An Accurate Numerical Method and Algorithm for Constructing Solutions of Chaotic Systems[END_REF][START_REF] Pchelintsev | On the Poisson Stability to Study a Fourth-Order Dynamical System with Quadratic Nonlinearities[END_REF], a software in the C++ language for Linux was developed. The source codes are available on GitHub [START_REF] Pchelintsev | The Reliable Calculations for the 4-th Order System[END_REF].

To perform matrix operations in C++, the algorithms and template types of the uBLAS library of the linear algebra collection of the Boost class libraries were used. Note that the uBLAS library gives a large gain in time.

Usually, in calculations, many researchers work with single or double precision subnormal real numbers represented in the IEEE 754 format. The main disadvantage here is the fixed accuracy of the representation of real numbers which may not allow us to numerically construct approximations to unstable solutions of systems of differential equations over long time intervals. Therefore, for high-precision calculations, the MPFR C++ library [START_REF]GNU MPFR Library for Multiple-Precision Floating Point Computations with Correct Rounding[END_REF] is used.

Algorithm 1 was applied to check the accuracy of the found unstable cycle [START_REF] Pchelintsev | A Numerical-Analytical Method for Constructing Periodic Solutions of the Lorenz System // Differencialnie Uravnenia i Protsesy Upravlenia[END_REF] in the system (2):

x 1 (0) ≈ -2.147367631, x 2 (0) ≈ 2.078048211, x 3 (0) = 27. ( 9
)
The period value is obtained equal to T c ≈ 1.558652210. The initial values [START_REF] Liao | On the Clean Numerical Simulation (CNS) of Chaotic Dynamic Systems[END_REF] were checked on the period in the computer program [START_REF] Pchelintsev | The Reliable Calculations for the 4-th Order System[END_REF] that implements the FGBFI method with ε pw = 10 -25 , 100 bits for mantissa of real number and ε m = 1.57772 • 10 -30 . With such parameters of the method, the approximate values x 1 (T c ), x 2 (T c ) and x 3 (T c ) were also verified by the same numerical method, but in backward time. The values at the backward time in the end of the trajectory arc coincide with (9) up to the 9th character inclusive after the point. The resulting values of x 1 (T c ), x 2 (T c ) and x 3 (T c ) coincide with (9) up to the 8th character inclusive.

The cycle corresponding to ( 9) is shown in Figure 1.

- D. Viswanath [START_REF] Viswanath | The Fractal Property of the Lorenz Attractor[END_REF] found in the system (2) an unstable cycle with a long period and with 99 decimal places use the Lindstedt-Poincaré technique:

x 1 (0) ≈ -13.568317317591138693791116532738086146665425413802770267307341928920639925115986035379124247913350182,

x 2 (0) ≈ -19.134575113926898610482106145675590555238063694018831440659257986585209042730623744601562225619287641, x 3 (0) = 27, T c ≈ 171.86372913973174676014481271369986804353836527546572814842984169003209300163561597123596376874805444. ( 10
)
The author of this article decided to check the found values with high accuracy in the forward time by the FGBFI method. For this, it was taken 390 bits for mantissa of real number. Then ε m = 7.93107 • 10 -118 . The value of accuracy for the power series is ε pw = 10 -110 . The maximum degree of the approximating polynomial at variable integration steps is 78. The computation time was approximately 6.2 min.

After the computational experiment, it was found that at the end of the period the values x 1 (T c ), x 2 (T c ) and x 3 (T c ) coincide up to 38 decimal places with [START_REF] Lozi | A New Reliable Numerical Method for Computing Chaotic Solutions of Dynamical Systems: the Chen Attractor Case[END_REF]. At the same time, the increase in accuracy (decreasing values ε m and ε pw ) does not affect the obtained results where we can conclude that some decimal places of the values [START_REF] Lozi | A New Reliable Numerical Method for Computing Chaotic Solutions of Dynamical Systems: the Chen Attractor Case[END_REF] are incorrect in [START_REF] Viswanath | The Fractal Property of the Lorenz Attractor[END_REF], or need more correct decimal places. Note that the classical numerical methods (such as the Euler method, Runge-Kutta methods, Adams methods, etc.) do not make it possible to carry out such a high-precision check of the values [START_REF] Lozi | A New Reliable Numerical Method for Computing Chaotic Solutions of Dynamical Systems: the Chen Attractor Case[END_REF] due to the instability of the found cycle by D. Viswanath.

We also note that the use of high-precision calculations is important in data encryption. For example, the results obtained in [START_REF] Flores-Vergara | Implementing a Chaotic Cryptosystem in a 64-bit Embedded System by Using Multiple-Precision Arithmetic[END_REF] show that when using the high-precision arithmetic, the generated sequences provide good randomness and security with more iterations of the implemented discrete-time chaotic systems compared to the received sequences in the IEEE 754 format.

The Systems of Explosive Type and Application of the FGBFI Method

The dynamic system of the explosive type is a system for which the state ∥X(t)∥ → ∞ is reached in a finite time. The examples of such systems are the system (3), and the population explosion model

ẋ1 = kx 2 1 ,
where the coefficient k > 0.

In other words, in this section of the article we will consider the autonomous systems of differential equations whose solutions have the vertical asymptotes. Such systems have been little-studied in the literature. We note the papers [START_REF] Zhukovskii | On a Parametric Specification of the Solution of a Differential Equation and Its Approximate Construction // Russian Mathematics[END_REF][START_REF] Zhukovskaya | About Vertical Asymptotes of Integral Curves of Ordinary Differential Equations[END_REF] where the scalar case is described: the asymptote existence theorem is proved and a modification of the Euler method is given which allows one to localize the asymptotes of solutions without jumping them.

For the multidimensional case, no studies have been found in the known literature. Most likely, this is due to the fact that for these systems, at ∥X(t)∥ → ∞, the energy also tends to infinity (which cannot be in real systems according to the energy conservation law); however, the systems of explosive type can be models of real systems in some approximation. Therefore, it is important for the researcher to have an idea at what time t as a sharp increase of phase coordinates will occur when it is not possible to find solutions of the system in quadratures.

Before proving the theorem on the existence of vertical asymptotes, we prove the following statement:

Lemma 1. Let the number sequence r (1) , r (2) , ... converge to some number r, the sequence w (1) , w (2) , ... converge to number w, and there is the number d > 0 that r (m)w (m) ≥ d

for all m ∈ N, then r > w.

Proof. Let us represent the sequences as

r (m) = r + ζ (m) , w (m) = w + ϑ (m) ,
where ζ (m) and ϑ (m) are the infinitesimal sequences.

From the inequality [START_REF] Lozi | A New Accurate Numerical Method of Approximation of Chaotic Solutions of Dynamical Model Equations with Quadratic Nonlinearities // Chaos[END_REF] it follows that

r (m) -w (m) = (r -w) + ζ (m) -ϑ (m) ≥ d > 0. ( 12
)
Note that the sequence

ζ (m) -ϑ (m)
is also infinitesimal. Therefore, there exists a number m = m * such that for all m > m *

ζ (m) -ϑ (m) ≤ ε,
moreover, we choose the small positive number ε such that

d -ε > 0, because the number d is fixed. Then from (12) (r -w) + ε ≥ d and r -w ≥ d -ε > 0,
whence it follows that r > w.

Note that the proof of this lemma can be carried over to the case when the sequences r (m) (t) and w (m) (t) are the sequences of scalar functions that converge uniformly over some time interval.

Let us now prove a theorem that allows us to establish the existence of the vertical asymptotes in the system (1) at t > 0.

Theorem 1. If some equation with the number q of the system (1) can be represented as (k is some number of the phase coordinate

x k , k ̸ = q) ẋq = c q + ⟨ Qq ( X -A q ), X -A q ⟩ + υx 2 k , (13)
where the symmetric matrix Qq is positive definite (D + Q ) or negative definite (D - Q ) and is obtained by deleting a row and a column with number k from matrix Q q ,

X(t) = [x 1 (t) . . . x k-1 (t) x k+1 (t) . . . x n (t)] T , X ∈ R n-1 ,
A q = a q,1 . . . a q,k-1 a q,k+1 . . . a q,n T , the numbers υ > 0, c q ≥ 0 for D + Q or the numbers υ < 0, c q ≤ 0 for D - Q , the initial conditions x k (0) ̸ = 0, x q (0) > a q,q for D + Q or x q (0) < a q,q for D - Q , then the q-th component of the solution of the system (1) has a vertical asymptote at t > 0.

Proof. For definiteness, we will consider the situation D + Q . A similar proof can be carried out for D - Q . Then the inequality ( [START_REF] Demidovich | Lectures on mathematical theory of stability[END_REF], pp. 34, 35)

⟨ Qq ( X -A q ), X -A q ⟩ ≥ λ q ∥ X -A q ∥ 2 2 is valid where λ q > 0 is the smallest eigenvalue of the matrix Qq . Since

∥ X -A q ∥ 2 2 = n ∑ p=1 p̸ =k (x p -a q,p ) 2 ≥ (x q -a q,q ) 2 , then c q + ⟨ Qq ( X -A q ), X -A q ⟩ ≥ c q + λ q (x q -a q,q ) 2 . ( 14
)
Let us introduce the equation ẏ = c q + λ q (ya q,q ) 2 (15)

with the initial condition y(0) for which x q (0)a q,q > y(0)a q,q > 0.

By the existence and uniqueness theorem, the solution of the Equations ( 13) and ( 15) exists and is unique on some time interval [0, t * ].

Let us construct the integral equations equivalent to the Equations ( 13) and ( 15):

x q (t) = x q (0) + t 0 c q + ⟨ Qq ( X(η) -A q ), X(η) -A q ⟩ + υx 2 k (η) dη, y(t) = y(0) + t 0 c q + λ q (y(η)a q,q ) 2 dη.

(

) 17
Since the function x k (t) is continuous on the interval [0, t * ], x k (0) ̸ = 0 and υ > 0, then

R(t) = υ t 0 x 2 k (η)dη > 0 for t ∈ (0, t * ]. Let R * = R(t *
) is a maximum value of the function R(t) on the segment [0, t * ].

Let us write the function x q (t) as

x q (t) = R(t) + x q (0) + t 0 c q + ⟨ Qq ( X(η) -A q ), X(η) -A q ⟩ dη.

Let us now consider the schemes of successive approximations obtained from the Formula (17):

x (m+1) q (t) = R(t) + x (m) q (0) + t 0 c q + ⟨ Qq ( X(m) (η) -A q ), X(m) (η) -A q ⟩ dη, y (m+1) (t) = y (m) (0) + t 0 c q + λ q (y (m) (η) -a q,q ) 2 dη,
where the vector function X(m) (η) is obtained from X(η) by replacing the function x q (η) by x

(m) q (η). Let x (0) q (t) ≡ x q (0), y (0) (t) ≡ y(0).
Since the existence and uniqueness theorem is proved by the method of successive approximations, then for the constructed approximations

lim m→∞ x (m) q (t) = x q (t), lim m→∞ y (m) (t) = y(t)
uniformly on the time interval [0, t * ].

By virtue of the inequality [START_REF] Pchelintsev | The Reliable Calculations for the 4-th Order System[END_REF], at the first iteration we have

x

q (t) ≥ R(t) + x q (0) + t 0 c q + λ q (x q (0)a q,q ) 2 dη.

From the positiveness of the function R(t), the inequalities ( 16) and

(x q (0)a q,q ) 2 > (y(0)a q,q ) 2 , we have

x

q (t) > y(0) + t 0 c q + λ q (y(0) -a q,q ) 2 dη ≡ y (1) (t), (1)
i.e., x

q (t) > y (1) (t) for t ∈ [0, t * ]. (1)
Also from the the inequalities ( 16) and ( 18), we have

x

q (t)a q,q > y (1) (t)a q,q > 0.

Let x (m) q (t)a q,q > y (m) (t)a q,q > 0 for t ∈ [0, t * ]. Hence we get x (m) q (t) > y (m) (t). In a similar way, one can show that x (m+1) q (t)a q,q > y (m+1) (t)a q,q > 0 and x (m+1) q (t) > y (m+1) (t).

At that max

t∈[0,t * ] x (m+1) q (t) -y (m+1) (t) ≥ R * > 0 for all m ∈ N.
Then, by virtue of the proved lemma x q (t) > y(t)

at t ∈ [0, t * ].

Note that the solution to the Cauchy problem ( 15)-( 16), which is easy to find in quadratures, has a vertical asymptote. Then it follows from the inequality (19) that the function x q (t) also has a vertical asymptote. Let us consider an example of a system that satisfies the conditions of this theorem. In the article [START_REF] Gotthans | [END_REF], the following system was introduced:

       ẋ1 = ax 3 , ẋ2 = x 3 f 1 (x 1 , x 2 , x 3 ), ẋ3 = x 2 1 + x 2 2 -r 2 + x 3 f 2 (x 1 , x 2 , x 3 ), (20)
in which the equilibrium positions lie on the circle x 2 1 + x 2 2 = r 2 in the plane x 3 = 0, a is the parameter of the dynamical system. Usually the functions f 1 and f 2 are chosen as

f 1 (x 1 , x 2 , x 3 ) = bx 1 + βx 2 3 , f 2 (x 1 , x 2 , x 3 ) = γx 1 ,
where b, β and γ are also system parameters.

Let r = 0, β = 0 and f 2 (x 1 , x 2 , x 3 ) = αx 3 , where α > 0 is a parameter. We get

       ẋ1 = ax 3 , ẋ2 = bx 1 x 3 , ẋ3 = x 2 1 + x 2 2 + αx 2 3 . (21)
The system [START_REF] Demidovich | Lectures on mathematical theory of stability[END_REF] belongs to the class of systems (1) and has the unique equilibrium (0, 0, 0). For the last equation of the system c 3 = 0, Q3 = 1 0 0 α , A 3 = 0, the matrix Q3 is symmetric and positive-definite. Then the component x 3 (t) of the solution of the system (21) has a vertical asymptote for x 3 (0) > 0 and x 1 (0) ̸ = 0. Note that the authors of the article [START_REF] Gotthans | [END_REF] do not provide a qualitative analysis of the system (20). As shown above, it may happen that it has unbounded solutions. Therefore, before looking numerically for attractors in dynamical systems, the researcher needs to show the existence of bounded solutions, for example, using the Lyapunov functions.

Note that by the proof of the above theorem, it is not necessary that all equations of the system (1) have a quadratic right-hand side. It suffices to have at least one equation that satisfies the conditions of this theorem.

It follows from the Formula (7) that the application of the FGBFI method allows one to approach the vertical asymptote arbitrarily close without jumping over it, since τ(Υ 0 ) = O 1 ∥Υ 0 ∥ 2 at ∥Υ 0 ∥ → ∞, i.e., the Formula (7) give a guaranteed region of convergence of the power series, and for larger ∥Υ 0 ∥ the integration step is smaller. Thus, from step to step applying Algorithm 1, we will come as arbitrarily close to time t as localizing the moment of explosion.

7 :

 7 ∆t := Tct 8: ct := T 9: else 10: ct := ct + ∆t 11: end if 12: Υ 0 := X l-1 13: X l := X l-1 14: p := 1 ▷ The product of powers of t 15: i := 0 ▷ Number of the current member of the series 16: do 17: i := i + 1 18:p := p • way • ∆t ▷ Calculate the current degree of t = ∆t ▷ taking into account the direction of integration; p is negative for odd powers ▷ when way = -119:

the current member of the series 21 :

 21 L := ∥Υ i ∥ • |p|▷ The right side of inequality (8) 22: while L > ε pw 23: Print the value of time way • ct and the vector X l of the received phase coordinates 24: if X l / ∈ S a then 25:

Figure 1 .

 1 Figure 1. The cycle in the system (2) obtained by the FGBFI method.

end if 32: if ended then 33: Terminate the algorithm 34: end if 35:

	27:	ended := true
	28: end if
	29: if ct = T then
	30:	ended := true
	31:

Note also that the application of the FGBFI method for the numerical integration of the equation [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF] gives an approximate value of the number π localizing the asymptote (4).

The Piecewise Smooth Systems with Quadratic-Type Nonlinearities. The Systems with Hysteresis

In practice, the problem often arises of studying piecewise smooth systems, for example, the Chua system [START_REF] Kuznetsov | Hidden Attractors Localization in Chua Circuit via the Describing Function Method[END_REF]:

where ρ, θ, χ, m 0 and m 1 are system parameters, the function

The authors of the article [START_REF] Kuznetsov | Hidden Attractors Localization in Chua Circuit via the Describing Function Method[END_REF] use the method of a describing function to localize hidden attractors in the system [START_REF] Gotthans | [END_REF], and the values of the initial conditions are also given. Therefore, the task of developing the high-precision numerical methods for constructing approximate solutions to systems of the form ( 22) is important. Note that the described idea can also be transferred to models with hysteresis which were considered in [START_REF] Medvedskii | Unstable Oscillating Systems with Hysteresis: Problems of Stabilization and Control[END_REF][START_REF] Meleshenko | Conservative Chaos in a Simple Oscillatory System with Non-Smooth Nonlinearity[END_REF][START_REF] Semenov | Self-Oscillations in a System with Hysteresis: the Small Parameter Approach[END_REF].

We can extend the FGBFI method to piecewise smooth systems with quadratic nonlinearities. Let us describe the idea of the method using the system (22) as an example.

From the form of the function ϕ(x 1 ), the phase space is divided into three parts, where the right side of the ( 22) system is smooth. Similarly, the phase space can be divided into regions of smoothness for systems with hysteresis. For example, for systems with a non-smooth potential [START_REF] Meleshenko | Conservative Chaos in a Simple Oscillatory System with Non-Smooth Nonlinearity[END_REF]. The joint point coordinates are

By the initial condition for the coordinate x 1 , we determine in which part of the phase space the initial point is located. The modules in the Formula (23) expand accordingly. Next, we take a forward step in time ∆t according by Algorithm 1. In this case, it is necessary to remember the obtained polynomials x1 (t), x2 (t) and x3 (t) approximating the corresponding phase coordinates x 1 (t), x 2 (t) and x 3 (t) on the time interval [0, ∆t]. If x1 (∆t) does not belong to the current part of the phase space, then it is necessary to find the time t pl with a high accuracy when the trajectory intersects one of the planes [START_REF] Medvedskii | Unstable Oscillating Systems with Hysteresis: Problems of Stabilization and Control[END_REF]. For this, for example, one of the equations x1 (t pl ) = ±1 is numerically solved by the secant method.

Further, the vector of initial conditions is taken equal to

and for the new part of the phase space, the modules are again expanded in the Formula ( 23). The advantage of the described algorithm is that we find the moments of intersection of the trajectory of the system [START_REF] Gotthans | [END_REF] with the planes [START_REF] Medvedskii | Unstable Oscillating Systems with Hysteresis: Problems of Stabilization and Control[END_REF] with a high accuracy, slightly accumulating calculation errors.

Conclusions

In this article, the high-precision FGBFI method for constructing approximations to unstable solutions of dynamic systems of the form (1) was described, and a link to the developed program software was also provided. The arbitraryprecision numbers were represented by mpreal data type of the MPFR C++ library with overloaded arithmetic operations, as well as friendly mathematical functions.

For systems of explosive type, a theorem is proved that allows one to establish the existence of asymptotes in the system (1), and the advantage of the numerical FGBFI method in this case is also shown.

The FGBFI method has also been extended to piecewise smooth systems with quadratic-type nonlinearities using the Chua system as an example, as well as systems with hysteresis.


Ramzi Chaieb
email: ramzi.chaieb@hotmail.com

Karimkalti
email: karim.kalti@gmail.com

Najoua Essoukri

Ben Amara

Interactive Content-Based Document Retrieval Using Fuzzy Attributed Relational Graph Matching

Keywords: Fuzzy Attributed Relational Graph indexing, Graph matching distance, Document image retrieval, User interaction

The evolution of digital technology and the desire to maintain an easy access to scanned documents implies the interest of Document Image Retrieval (DIR). In this paper, we propose a fuzzy graph-based document retrieval approach for determining similarity between document images for content based document retrieval. To model the structure of document images, we have opted to use Attributed Relational Graphs (ARGs). For each document region (text, graphic, etc) we associate a node which is characterized by a set of fuzzy membership degrees reflecting low-level properties (texture, shape, color, etc). We use this fuzzy description in order to guarantee more robustness against the eventual segmentation errors which may be occurred after the segmentation of document regions. ARGs edges represent spatial relationships between regions which have common boundaries. Finally, we have developed a tree-search based optimal matching algorithm, which allows the search for document according to its structure.

The database used for experiments is composed of segmented Coran document images from the National Library of Tunisia.

Different weights are assigned for regions and edges according to their relative importance. The results obtained demonstrate the effectiveness of fuzzy graph-based document retrieval. Such an approach is very useful for several applications in many fields.

I. INTRODUCTION

In recent years, with the rapid development of digital devices for image creation, storage and transmission, huge amounts of images are produced every day (for example in administrations, museums, etc). The huge amount of scanned documents stored on the internet is growing rapidly. Manual organization of documents is too costly and sometimes even impossible. Automated content-based document management methods, generally known as DIR techniques, are needed to deal with these amounts of data.

Motivated by the enormous demand of information accessing in document image databases, more and more research efforts have been devoted to DIR in recent years [1, [START_REF] Pattanaik | Beginners to Content Based [mage Retrieval[END_REF][START_REF] Deb | An Overview of Content-based Image Retrieval Techniques[END_REF]. In DIR systems, the user specifies a query either in the form of image, text or sketch and the system is expected to return the items that are visually and semantically most similar to the query. Traditional approaches for D1R are text-based [START_REF] Pattanaik | Beginners to Content Based [mage Retrieval[END_REF][START_REF] Sameriya | A Survey on Content Based Image Retrieval (CBIR) Schemes[END_REF][START_REF] Ahmed | A Study on Different Image Retrieval Techniques in Image Processing[END_REF][START_REF] Datta | Image Retrieval: Ideas, Influences, and Trends of the New Age[END_REF]. A key task of text-based approaches is to annotate and 978-1-4799-1805-8/15/$31.00 ©2015 IEEE index the images with textual descriptions or keywords for retrieving images [START_REF] Pattanaik | Beginners to Content Based [mage Retrieval[END_REF][START_REF] Sameriya | A Survey on Content Based Image Retrieval (CBIR) Schemes[END_REF][START_REF] Ahmed | A Study on Different Image Retrieval Techniques in Image Processing[END_REF][START_REF] Datta | Image Retrieval: Ideas, Influences, and Trends of the New Age[END_REF]. However, many limitations of the text-based approaches make them still far from working in real-world applications. Describing the rich content of a document image by only a few keywords is almost difficult. In addition, text descriptions reflect the subjectivity of the annotator and the annotation process is prone to be incomplete, ambiguous and very difficult to be automated. Hence, fully text-based approach is not practical enough for DIR applications. This has created a demand for effective and flexible techniques for automatic DIR to overcome the shortcomings of the text-based retrieval mechanism. Content Based Document Retrieval (CBDR) has been suggested as an alternative approach for information retrieval to find document images which are visually most similar to a given query from a document image database. At present, CBDR has become an extremely active area of research in the visual information retrieval [START_REF] Zhao | Content-Based Image Retrieval Using Optimal Feature Combination and Relevance Feedback[END_REF][START_REF] Singhai | A Survey On: Content Based Image Retrieval Systems[END_REF]. In CBDR, the document images in the database are indexed by extracted visual features (such as color, shape, texture, etc.). Based on the indexing scheme with visual descriptions of low-level features, visual queries can be formulated and similarity of document images can be measured by employing a distance function defined in the CBDR systems [1,[START_REF] Pattanaik | Beginners to Content Based [mage Retrieval[END_REF][START_REF] Deb | An Overview of Content-based Image Retrieval Techniques[END_REF].

A document image is generally composed of several components such as logo, text, figure, table, signature, etc. Usually, document images belonging to the same type have similar layout structures. For example, a letter is generally composed of a letter head, text area and a signature. An article contains columns of text blocks, figures, tables, etc. Traditional CBDR approaches are not powerful enough to capture these high-level concepts because of the semantic gap between high level concepts and low-level features. A segmentation phase is usually necessary to separate adjacent regions. However, segmentation errors may occur frequently. From the extracted regions, a structural signature is needed to define document image regions and their topological relationships. Graph based approaches usually provide a rich and holistic description of the layout and content of the analyzed document images. In the literature, a lot of research has been done on graph-based representation and graph matching algorithms [START_REF] Sameriya | A Survey on Content Based Image Retrieval (CBIR) Schemes[END_REF][START_REF] Datta | Image Retrieval: Ideas, Influences, and Trends of the New Age[END_REF][START_REF] Sharma | Determining similarity in histological images using graph-theoretic description and matching methods for content-based image retrieval in medical diagnostics[END_REF][START_REF] Philipp-Foliguet | FReBIR: An image retrieval system based on fuzzy region matching[END_REF].

In this study, we propose a fuzzy graph-based document image representation in order to reduce the problem of possible segmentation errors. Each segmented region is characterized by a set of fuzzy membership degrees derived from low-level features (texture, shape, color, etc). For example, a segmented region may be a 80% text and 20% graphic. In addition, a similar set of fuzzy membership degrees is used to specify the spatial relationship (above, below, right, left, surrounded-by, etc) between two adjacent segmented regions. Another contribution of this work is the human's interaction in order to surmount the limited capacity of traditional CBDR systems. User's weights are involved in the retrieval process.

The remaining of this paper is organized as follows. The next section gives an overview of the proposed approach. Section 3 presents the document F ARG representation. Sections 4 describes the process of F ARG matching algorithm. Experimental results are evaluated and discussed in section 5. Finally, concluding remarks are given in the last section.

II. METHODOLOGY

The proposed CBDR approach is composed of two main phases: indexing and retrieval phases (Fig. 1).

Document Image

A. Offline Indexing Phase

The used document images are indexed into a database. Each document image is represented by a F ARG that describes the document image as a set of regions with relationships between them. Both the nodes and the edges are represented with their attributes.

B. Online Retrieval Phase

Once a user specifies an example document image, the proposed query is then represented by a FARG. The retrieval task is to fmd the closest matches of the query F ARG representing the structure of a document image components to a database of F ARGs already generated. A similarity measure is used in the comparison of two document images. The obtained costs of graph matching are used to sort the database document images in a decreasing order of similarity. Finally, the obtained results are ranked as the fmal retrieval result and returned to the user.

III. Fuzzy ATTRIBUTED RELATIONAL GRAPH REPRESENTATION

We have used FARGs as a structural representation of document layout. One motivation for the F ARG representation is to reduce the effect of inaccurate segmentation on retrieval.

F ARGs have an advantage over vector-based features as they can model the spatial relationship between different regions [START_REF] Sharma | Determining similarity in histological images using graph-theoretic description and matching methods for content-based image retrieval in medical diagnostics[END_REF][START_REF] Philipp-Foliguet | FReBIR: An image retrieval system based on fuzzy region matching[END_REF].

A. Attributed Relational Graph

An ARG is a graph whose nodes and edges are both represented with attributes. It can be represented as a 3-tuple (N; E; A) where N is the set of nodes (vertices), E � N x N is the set of edges (relations) and A is the set of attributes. Each node (n) attributes are defined as a vector Xi=[X/I)j, (u= 1,2,3'00" U), where U is the total number of node attributes associated with node ni' Similarly, for each edge (ej) a set of attributes are represented as a vector Yj=[Yj(v)j, (v=1,2,3,00., V), where V is the total number of edge attributes associated with edge ej. An example of an ARG is shown in Fig. 2.

Xl � [X//). • .... X/!)} ..... Y/'? X2 � [X//). .. .. X/'l } X3 � [x.tI, .
. , xl!)} Fig. 2. An example of an ARG with node attribute vectors {Xi, i � I, 2, 3} and edge attribute vectors (Y J ' j � 1, 2, 3) B. Regional Adjacency Graph Regional Adjacency Graph (RAG) is an ARG whose nodes represent document image regions and edges represent spatial relationships between regions which have common boundaries. Each node has a unique identifier number. Node attributes represent informations about the region such as label, type of the region (which may be text, graphic, etc), area, perimeter, shape, color, etc. Edge attributes specify relationships between two regions like the distance between their centroids, common boundary length, orientation which may be above, below, right, left, surrounded-by, etc. An example of a RAG is shown in Fig. 3. One advantage of F ARG is the representation of document images with a fuzzy approach, which is similar to human perception [10, ll]. Fuzzy logic approaches are proposed to map from the low-level features to high-level concepts in order to bridge the semantic gap between them. Document image regions often possess properties that are fuzzy in nature and it has been extended to include fuzzy information into the attributes.

[n a FARG, each vertex may take attributes from the set Z={z; I i = t, 2, ... , f}. For each attribute Z;, it will take values from S; = {sij I j = I, 2, ... , J;}. The set of all possible fuzzy attribute-value pairs is Lv = {(z; , AsJ I i = 1, ... , fj, where As; is a fuzzy set on the attribute-value set S;. A valid pattern primitive is just a subset of Lv in which each attribute appears only once. The set of all those valid pattern primitives is denoted as II. Each vertex will be represented by an element of fl.

Similarly, each arc may take attributes from the set F={f;Ii=l, 2, ... , 1} in which each j; may take values from T; = {tij I j = 1, 2, ... , .!;} and La = {(f; , BTJ I i = 1 , ... , f} denotes the set of all possible relational attribute value pairs, where Bn is a fuzzy set on the relational attribute value set T;. A valid relation is just a subset of La in which each attribute appears only once. The set of all those valid relations is denoted as CPo A Fuzzy Attributed Graph (FAG) over L=(L"LcJ, with an underlying graph structure H=(N,£) is defined to be an ordered pair (V,A), where V= (N,(j) is called a fuzzy vertex set and A =(£,6) is called a fuzzy arc set. The mapping (j: N -7[1 and 6:£ -7<P are called fuzzy vertex interpreter and fuzzy arc interpreter, respectively. This definition also applies when there are non-fuzzy attributes, since a crisp (non-fuzzy) set can always be represented as a special case of a fuzzy set.

IV. Fuzzy ATTRIBUTED RELATIONAL GRAPH MATCHING ALGORITHM

[n this section, a fuzzy graph matching algorithm is used to calculate the graph matching distance between a query document image and a document image in the database. Given a query document image, its F ARG is matched to each F ARG in the database. A similarity measure between matched FARGs is assigned to each pair of FARGs.

Let G(N,£,a,/3) and G'(N',£',a',/3') be respectively a query F ARG and a target F ARG. For each node n; (i E 1, 2, ... , N) in G, the set of nodes membership degrees is given as a;: a ; =a ; ( 1 ),a/2), ... ,a/a),a E {1, 2, ... ,A} (1

)
where N is the number of nodes in G and A is the total number of nodes membership degrees associated with each node n;. We denote by a the set of all vectors a;.

(2)

Nodes membership degrees indicate how much a region belongs to a graphic zone, text zone, etc. Different weights are assigned for each node features according to their relative importance to describe regions. If we denote by wn i the weight of node n;, the query image can be represented as a weight vector W N in the weight space by: N WN = [wnl, wnz' ... , wn N ], i E {1,2, ... , N} and L wn i = 1(3)

;= 1

Node n; with the highest value of WN is considered the most important node during the matching of nodes.

Similarly, for each edge e j 0 E I, 2, ... , £) in G, the set of edges membership degrees is given as A /3j = /3/ 1 ),/3/2), ... ,/3/b),b E {1, 2, ... ,B} [START_REF] Zhao | Content-Based Image Retrieval Using Optimal Feature Combination and Relevance Feedback[END_REF] where £ is the number of edges in G and B is the total number of edges membership degrees associated with each edge ej. We denote by /3 the set of all vectors /3j.

(5)

Edges membership degrees specifies the spatial relationship between two nodes (above, below, right, left, surrounded-by, etc). [f we denote by wei the weight of edge e;, the query image can be represented as a weight vector Wli in the weight space by: E WE = [Wei' Wez' ... , We E ]' i E {1,2, ... , E} and L We i = 1 (6) ;= 1 Edge e; with the highest value of W", is considered the most important edge during the matching of edges.

The graph matching is implemented through a tree-based search. Assume that G an G' with order m and n (m :s n) , respectively. [n the special case where m=n, the problem is to find the optimal isomorphism between G and G'. First, the decision tree is constructed. [t has height m, and n-p sons for each node at levelp = 0,1, ... , m. Therefore, at any levelp>O, there is a path that consists of ordinal nodes in the decision tree from the root to node N. The first level of the tree contains all possible combinations of matched nodes. A matching between a query node and a target node is permissible only if they have the same nature (text, graphic, etc). Then, the cost of each node pair is computed using weighted Euclidean distance (Equation 7).

Tn dn(qn,sn) = L Wn i X (qn i -Sn) 2 [START_REF] Ahmed | A Study on Different Image Retrieval Techniques in Image Processing[END_REF] where qn is a node of G, Sn is a node of G', Tn is the number of node membership degrees, qn i is a membership degree of the query node, sn i is a membership degree of the mapped node.

[n order to increase the accuracy of our approach, we introduced weighting wn i to our similarity measurement (wn i is the importance weight for a node n;). Node pair (q;,sJ with the lowest cost is expanded.

(8)

The cost of each edge pair is computed using weighted Euclidean distance (Equation 9). ; where qe is an edge of G, Se is an edge of G', Te is the number of edge membership degrees, qe i is a membership degree of the query edge, se i is a membership degree of the mapped edge. Similarly to node weighting, we introduced weighting we i to our similarity measurement (we i is the importance weight for an edge e;).

The tree is extended until path with least total cost giving optimal mapping is found. An example of the tree-based search approach is given in Fig. 4. The total node and edge distances are computed using equations 10 and 11, respectively.

N w TND = L dn i (qn,sn) (10) TED = L deJqe' se) (11) i=l i=l

where Nand Ware the number of matched node pairs and the number of matched edge pairs, respectively.

The total similarity measure between a F ARG pair can be defined by the formula 12:

SIM(G, G') = w G x TND + w G ' x TED ( 12
)
where W G and w G ' are properly selected weights of G and G'.

V. EXPERIMENTAL RESULTS

Experiments were conducted on a database composed of 200 initially segmented Coran document images from the National Library of Tunisia In the retrieval phase, firstly user assigns an importance degree to each region by giving an appropriate value between 0 and 1 for each region. Region with the highest value of importance degree is considered the most important region during the matching of regions. Similarly, user assigns an importance degree to each relationship. So, the search will rather focuses the most important relationships selected by the user. Assigning importance degrees provides a more precise formulation of user's queries. Figure 7 shows an example of the used importance degrees related to the nodes and edges of the query graph G. In this example, region I is considered the most important region in the query document image. Each image is represented by a F ARG. Then, images whose F ARGs are closest in distance to a given query F ARG are selected and the images belonging to each selected F ARG are retrieved and displayed to user.

In these experiments, each region i is described by two membership degrees aF) and ap) which corresponds respectively to the percentage of each region to be a graphic or a text. Similarly, each relationship j between adjacent regions is determined using two metrics f3r and f3F) which corresponds respectively to the distance and the horizontal angle between the two centroids of adjacent regions.

VI. CONCLUSION

In this paper, we have proposed a CBDR approach based on F ARG representation to address the problem of matching between document images. F ARG representation has helped in reducing the effect of inaccurate segmentation on retrieval. The proposed approach is composed of two main modules that work in pipeline mode, The first module consists of query construction and F ARGs generation. The second module consists of F ARG matching approach, A tree-based approach has been developed to calculate the matching cost between two FARGs. Finally, best matching results are retrieved and displayed to the user. The experimental results are promising. They can be improved, as future work, using relevance feedback techniques in order to learn from the user's evaluation.

Besides, the proposed approach can be generalized to many types of document images such as administrative documents, journals, books, etc.

Fig. l .

 l Fig. l. Block diagram of the proposed CBDR approach

Fig. 3 .

 3 Fig. 3. An example of a RAG generation, regions are represented by nodes C. Fuzzy Attributed Relational Graph

  Te de(qe,se) = L we i X (qe i -se Y[START_REF] Sharma | Determining similarity in histological images using graph-theoretic description and matching methods for content-based image retrieval in medical diagnostics[END_REF]

Fig. 4 .

 4 Fig. 4. An example of tree-based matching process

  [12]. This database is created in collaboration with the National Library of Tunisia [12]. The segmentation method is based on the exploitation of connected components and morphological operators. It allows the extraction of three types of regions (graphic, text and frame). The database is divided into 5 classes. Each class contains 40 Coran document images with similar structures and contents. Variability between classes is caused by differences between the regions' specific characteristics and their positions in the Coran document images. Five examples from the five different classes are shown in Fig. 5. Two sample images from the used database and their segmented regions are shown in Fig. 6.

Fig. 5 .

 5 Fig. 5. Five examples from the five used classes [12]: (a) containing only text regions, (b) containing text regions surrounded by a frame, (c) containing graphic regions followed by text regions, (d) containing text, graphic and frame regions, (e) containing only graphic regions

Fig. 6 .

 6 Fig. 6. Two sample images from the database used in experiments [12] As shown in Fig. 7, Document image 1 and Document image 2 are represented by two graphs G and G: respectively. Their tree-based matching process is illustrated in the same figure. The total similarity measure between the two graphs G and G' is computed as follows: TND = dn5 (of Level 0) + dn5 (of Levell) + dnl (of Level 2) = 0.0077+ 0.0219 + 0.0379 = 0.0675 TED = del + de2 = l.8828 + l.0678 = 2.9506 SIM(G,G') = We; X TND + wG'x TED = 0.7 X 0.0675 + 0.3 X 2.9506= 0.9324 Precision and recall are the most popular measures used for evaluation of Content-Based Image Retrieval system performance [8, 9]. Precision and Recall are calculated only for one query. To evaluate a system for several queries, we use the mean precision and the mean recall. Examples of recorded Precision and Recall values are given in table 1. The obtained results are satisfactory. They are related to the number and the nature of regions (text, graphic or frame) and their relative positions in each Coran document image.

Fig. 7 .

 7 Fig. 7. An example of two FARGs (G and G') matching which corresponds to "Document image I" and "Document image 2", respectively

TABLE L OBTAINED

 L RESULTS (PRECISION, RECALL AND MATCHING TiME) REFERENCES

	Query Number	Relevant Retrieved Images	Retrieved Images	Relevant Images	Precision (%)	Recall (%)	Matching Time (s)
	I	26	34	40	78.47	65	0.0108
	2	25	32	40	78.12	62. 5	0.0134
	3	29	38	40	76.31	72. 5	0.0174
	4	26	35	40	74.28	65	0.0117
	5	27	37	40	72.97	67. 5	0.0142
	6	27	38	40	71.05	67. 5	0.0157
	7	29	41	40	70.73	72. 5	0.0153
	8	31	44	40	70.45	77. 5	0.0143
	9	33	48	40	68.75	82. 5	0.0162
	10	34	51	40	66. 66	85	0.0118
	Average				72.77	71.75	0.0141


Lionel Brunie

Laurent Lef

Laurent Lefevre

A DSM-based structural programming environnment for distributed and parallel processing

The paper describes an original programming environnment based on the DOSMOS 1 system.

Its implementation is built on a structural approach of parallel programming. By combining this structural model with weak consistencies schemes and protocols, we i m p r o ve the performances of DSM system and provide an original programming model that mixes message-passing and shared distributed objects. Morevoer by i n tegrating an easy development platform and montoring facilities, the DOSMOS system has been designed to provide a performant user-friendly programming environment.

Introduction

If distributed memory MIMD machines allow v ery high performances, their programming remains "esoteric" for most of end-users accustomed to classical mono-processor programming. Oppositely, shared memory parallel computers are easier to program but badly adapted to applications generating many memory accesses. In this context, the purpose of Distributed Shared Memory systems (DSM) is to implement, above a distributed memory architecture, a programming model allowing a transparent manipulation of virtually shared data. Thus, in practice, a DSM system has to handle all the communications and to maintain the coherence of the shared data.

In that framework, this paper describes an original programming environment, called DOSMOS 1 system. This system is based on a structural approach of parallel programming. In other words, DOSMOS proposes to the user to hierarchically structure the processes into groups and sub-groups of processes sharing a same set of variables. This feature, combined with optimized weak consistency protocols allows to reduce the amount o f c o m m unications required for the management of the shared data, and, as a consequence, to ensure e ciency and scalability to the applications. However, it would be unrealistic to argue that DOSMOS, or any other DSM system, e ciently deals with any kind of applications. That is why DOSMOS allows the programmer to mix both message-passing (PVM) and DOSMOS code. To complete the programming environment, DOS-MOS integrates a devoted monitoring tool (called DOSMOS-Trace) which has been added to the system to allow the user to understand the behavior of his applications.

At last, this programming environment has been designed to run both on distributed systems and on parallel machines. Thus, to ensure the portability of both the system and the applications, DOSMOS (as well as DOSMOS-Trace) has been developed on top of PVM.

1 DOSMOS : Distributed Objects Shared MemOry System 2 This paper is divided into three parts. After a short description of previous works (section 2), we analyse the basics of the DOSMOS DSM system in terms of management of shared data and process structuring (section 3). Then a description of the programming environment is proposed (development platform, programming model and monitoring facilities). At last, section 5 proposes a discussion both on the basic features of this programming environment and on implementation choices and points out future developments.

2 Purpose of Distributed Shared Memory systems and previous works By allowing the programmer to share "memory objects" (i.e. programming variables) in a transparent w ay, Distributed Shared Memory Systems (DSM) propose a interesting trade-o between the easy-programming of shared memory machines and the e ciency and scalability of distributed memory systems. Basically, a Distributed Shared Memory system is a mechanism that allows application processes to access to shared data in a transparent w ay. In other words, a DSM system releases the programmer from the management of all inter-process communications.

Both hardware and software implementations have been proposed. The main systems require to implement an additional software layer :

Virtual Shared Memory systems (VSM) allow to share pages of data, i.e. to merge into a single wide address space a set of memory pages distributed in the network. Such systems like MIRAGE FP89] or MUNIN CBZ91] h a ve to deal with speci c problems of operating systems.

Object-based Distributed Shared Memory systems (DSM) work at the program level, i.e. they implement a software layer that automatically generates, on the user's behalf, all the communications required to manipulate shared data. In other words, instead of de ning (and writing in the code) the inter-process communications, the programmer only speci es which data are actually shared. Then he can use these data as if they were local. On its side,the DSM system takes into charge all the communications necessary (as a message-passing programmer would do). Such DSM systems like O R CA TKB92] or CLOUDS system RAK89] have been implemented on parallel machines. The DOSMOS BL94, BL96] system belongs to this class of systems.

Basics of DOSMOS

DOSMOS is an object-based DSM system (cf section 2), i.e. it allows processes to share in a transparent w ay a set of passive objects (i.e. of programming variables) distributed in the network.

However, DOSMOS integrates novel features : DOSMOS Processes : Basically, a DOSMOS application is composed of two t ypes of processes: Application processes (A.P.) contain and execute the code (written in C) of the application Memory processes (M.P.) manage the whole DSM system, i.e. they provide A.P.

with the objects they request and maintain the data coherence. Each A . P . is connected to one and only one memory process. On the contrary, a n M . P . can be connected to several A.P. and several M.P.

Array allocation : DOSMOS allows to manipulate both basic type variables (integer, oat, char: : : ) and distributed arrays. These arrays are split into several \system objects", distributed in the network. Various splittings are provided : by r o w, by column, by block a n d by cyclic block. The system ensures a transparent access to arrays, whatever the splitting implemented.

Optimized weak consistency protocols : for e ciency and scalability purpose, DOSMOS allows to duplicate shared objects. It is clear that these replicas have t o b e k ept coherent.

Most of actually implemented models are strong consistency oriented. DOSMOS implements a w eak protocol: the release consistency. This model GLL + 90] provides two s y n c hronization operators: acquire and release. These operators allow processes which w ant to modify shared objects to lock and unlock them (in other words, these routines actually implement a m utual exclusion on the accesses to the shared objects).

Hierarchical structuring of the application processes : Previous DSM systems have a l w ays proposed \ at" models in which a n y shared object is accessible from any process. Such \anarchical" models cannot be scalable. In DOSMOS, processes can be group into groups in sub-groups in order to optimize the management of the coherence of data.

When one observes the behaviour of a DSM application, and more particularly the behaviour of a process participating to the application, it appears that if some shared data are intensively accessed by this process, some others are either very not often accessed or never accessed. This leads us to introduce some de nitions (see g. Usually, in previous systems, when an object O is modi ed, an invalidation message is sent t o all the processes P such that O 2 GV S (P). This prevents, as noted before, to ensure a good scalability. By using a hierarchical grouping of processes, DOSMOS limits the invalidation messages to processes such that O 2 LV S (P). 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 Basically, DOSMOS proposes to structure the application into hierarchical groups of processes sharing the same objects (Figure 2). In practise, a group is de ned by a set of processes and a set of shared objects. Processes of a same group share all the objects attached to the group, i.e. if they request an object, they will receive a copy of this object which will be automatically updated by the system.

11
But DOSMOS also allows processes to access to extra-group shared objects. For this purpose, in each group, a dedicated memory process, called Link Process (LP), plays the role of link between groups (see Fig 3). Thus, these special MPs takes into charge all the communications between groups. This model presents two important a d v antages :

the access to the shared objects is optimized the maintaining of the consistence is kept cheap. Experiments done with DOSMOS system on network of workstations and parallel machines have s h o wn great improvements when using hierarchical groups.

4 Programming environment

Development platform

The implementation of DOSMOS is based on two di erent l a yers (see Figure 4):

Preprocessing level: this layer analyses the user's application in order to detect and generate accesses to shared objects ( g. 5). This layer allows the system to be \transparent" DSM level: this layer assumes the creation and management of shared objects, groups and of various processes involved in execution of the application (APs, MPs, LPs). This management is performed using PVM routines.

in all the situations. Consequently, in order to allow the user to optimize speci c applications, DOSMOS allows to combine di erent programming models for user's confort. Consequently, three programming models are available:

Local programming : in order to minimize the accesses to shared objects, it is sometimes more performant t o w ork on local variables before modifying shared variables. D.S.M. programming : the user can use DOSMOS primitives to declare and access to the shared variables in a transparent w ay.

Mixing of DSM and message-passing programming : the user can integrate message-passing communications into DOSMOS applications. This feature presents two a d v antages First, it permits to deal with speci c applications. Second, it allows to port PVM applications on DOSMOS with slight modi cations of the code.

DOSMOS-Trace monitoring environment

The only way a user can in uence the behaviour of his application is the modi cation of the structure of the shared variables space. So, from the user point of view, monitoring facilities should allow him to precisely know the \activity" of the shared variables.

The purpose of the DOSMOS-Trace BLR96] monitoring environment is to provide such information in a scalable and weakly intrusive w ay. The DOSMOS-Trace tool is based on a set of dedicated processes which collect informations during execution. This data collection is completely transparent for the user. This tool provides several visualizations and informations about the execution like statistics on shared objects, histories: : : Such diagrams are extremely useful for the user to analyse problematical situations. Indeed they allow t o v ery easily isolate ping-pong e ects (e.g. g. 8), over-accessed variables, bottlenecks, not actually shared variables, etc. 5 Discussion and future works This paper has described a novel DSM-based programming environment, the DOSMOS system. In comparison with previous works, this system integrates original functionalities : structuring of the application processes into hierarchical groups, possibility of mixing message-passing code and DSM code, optimized weak consistency protocols, monitoring facilities.

The whole system has been designed to be as e cient and scalable as possible. Thus the process grouping allows, in conjunction with weak consistency protocols, to reduce the amount of communications required by the management of the DSM system. For the same purpose, the DOSMOS-Trace environment is based on distributed processes and has been designed to permit to use distributed traces les. Tests have s h o wn the e ectiveness of the approach developed in DOSMOS.

Opened to various programming models, designed to be e cient both on parallel machines and distributed systems, DOSMOS provides a portable development platform. Moreover by only adding few new primitives and by p r o viding graphical interfaces to analyse execution application, DOMOS is an user-friendly programming environment which can easily adapt to a non-expert parallel programming user.

We currently continue to improve our programming environment b y adding a new distributed tool to DOSMOS which will allow to debug the application code of the user. Moreover important DSM-applications (scienti c computing, imagery, neural networks: : : ) are currently being implemented on DOSMOS system to show the performance and e ectiveness of our approach.

  1) : G.V.S. : The Global Virtual Space (GVS) of a process is the set of the shared objects accessed (in read or write mod e ) b y a p r ocess during the execution of the application. L.V.S. : The Local Virtual Space (LVS) of a process is the set of the shared objects intensively accessed by this latter. E.V.S. : The Extern Virtual Space (EVS) of a process is the set of the shared objects rarely accessed b y t h e p r ocess. L e t P a p r o c e s s . W e h a ve : Global V irtual S pace(P ) = Local V irtual S pace(P ) + E xtern V irtual S pace(P )

Figure 1 :Figure 2 :

 12 Figure 1: Example of accesses distribution

Figure 3 :

 3 Figure 3: Groups and link processes

  Objet systeme: MAX(3) -Nombre et repartition des lectures

Figure 7 :Figure 8 :

 78 Figure 7: Number and origin of the read accesses performed on an object vs execution time (in black : i n ter-group accesses)

Programming model

As soon as the Use Dosmos() primitive has been executed, the user can access to the shared objects in a transparent w ay. H o wever, DOSMOS, as any DSM system, does not pretend to be e cient


S A Amin
email: s.a.amin@lut.ac.uk

D J Evans
email: d.j.evans@lut.ac.uk

L Lefèvre

PVM implementation for low-level image processing systolic array designs

This paper describes PVM implementation of systolic array designs of parallel algorithms for low-level digital image processing, in particular we consider the mean, weighted mean and laplacian filters. The PVM implementation of the design, is discussed with improvements of the original systolic design (macro and multi-pipeline). Comments and conclusions related to the implementation of the systolic array on PVM are provided in the performance section.

Introduction

The purpose of this paper is to identify a set of designs suitable for implementing low level image processing algorithms on distributed systems using Parallel Virtual Machine (PVM) model. We consider techniques for filtering digital images. Image smoothing has the purpose of removing noise and making the pixels' grey levels uniform. Most of the image processing algorithms need massive amounts of banded matrix operations. However, these algorithms contain explicit parallelism which can be efficiently exploited by processor arrays. All sections of the image have to be processed in exactly the same way, regardless of the position of the image section within the image, or the value of the pixel data. Low level functions involve matrix vector operations, which are repeated at very high speed. Typically, images must be processed in real time, at 25 images per second. Therefore with image sizes of 128x128, 256x256, and 512x512 or greater, there is a large amount of data to be processed in a highly repetitive process. Many low-level image processing operations, termed local operators, require access to the four or eight neighbouring intensity values of a pixel, when computing the new value for the pixel. Each member element in the image is replaced by some function of itself and the neighbouring elements within a window centred on that element. There are different architectures for carrying out these varied operations.

A number of systolic designs for 2D convolution and digital filters have been implemented [START_REF] Kung | A systolic 2d convolution chip[END_REF], [START_REF] Megson | An introduction to systolic algorithm design[END_REF], [START_REF] Wang | Gradient inverse weighted smoothing scheme and the evaluation of its performance[END_REF]. In this paper, a class of local smoothing filter algorithms are covered, mean, weighted mean and laplacian filters. For this purpose, we have developed several models of systolic arrays. Different kinds of cells are used. Various modifications of the systolic design are analysed, to handle a set of digital image filters.

The report starts with a brief definition of the algorithms. Another section presents a systolic array design for these filters.This paper also presents the implementation of a variety of digital image filter algorithms on the distributed computing using Parallel Virtual Machine (PVM) model [GBD + 93]. One of the aims was to design and build a programming workbench for developing image processing operations for low-level vision. The motivation for the work is to develop a methodology for the implementation of an image processing library on the distributed computing 2 Operator Algorithms.

We have developped various kinds of low-level image processing methods (like convolution, gradient filter, Inverse Gradient filter...). As example of our work, we present our systolic design of Laplacian Operator algorithms. The Laplacian operator is computed by convolving a mask with the image. The Laplacian low pass filter is a typical smoothing filter, a Laplacian mask is passed over the entire image, and the convolution operation is performed on each pixel. Each pixel is replaced by the sum of the products of the mask weighting and the appropriate neighbouring pixel values. Different choices are available when using this mask in two dimensions. The Laplacian filter using 3 by 3 convolution mask may be used. This utilises a mask or weighting matrix defined on two standard masks as shown in figure 1.

0 -1 0 -1 -1 -1 Laplacian 2 : -1 4 -1 Laplacian 3 : -1 8 -1 0 -1 0 -1 -1 -1 Figure 1: Laplacian masks
The value in the weighting matrix allowed a simpler and faster version of the algorithm than was obtained using the general convolution case. The algorithm given here 'smooths' a grey level input image and has I as an input image and O as an output image; both I and O contain M by M pixels, with P = M2. Each point of I for generating y, is as follows (for an image n by n ).

for i j = 1 2 ::::: n y ij = w 1 x i;1 j + w 2 x i j;1 + w 3 x i j + w 4 x i j +1 + w 5 x i+1 j By applying the weighting mask shown in figure 1 to this equation, then,

y i j = f (1 5) x i;1 j + x i j ;1 + x i j + x i j +1 + x i+1 j ]
Also if we employ the square window shown in figure 1 the output pixel is obtained by using the relation, y ij = w 1 x i;1 j;1 + w 2 x i;1 j + w 3 x i;1 j +1 + w 4 x i j;1 + w 5 x i j w 6 x i j +1 + w 7 x i+1 j;1 + w 8 x i+1 j + w 9 x i+1 j+1 and, y i j = f (1 9) x i;1 j ;1 + x i;1 j + x i;1 j +1 + x i j ;1 + x i j +x i j+1 + x i+1 j;1 + x i+1 j + x i+1 j+1 ]

This algorithm approach for the weighted mean filter is similar to the mean filter algorithm described in the previous section . The difference in this case, a weighted mean filter is often used in which the weight for a pixel is related to its distance from the centre point. For a 3x3 window, the filter weights are shown in figure 2. The neighbours that lie on the same side of the point (i,j) are weighted more heavily than the others.

1/6 1/16 1/8 1/16 Weighted mean 1 : 1/6 1/3 1/6 Weighted mean 2 : 1/8 1/4 1/8 1/6 1/16 1/8 1/16

Figure 2: Weighted mean masks

Systolic Array Design

The systolic array design to accommodate the previous filters, consists of k2 cells, where k2 is the size of the window. Each cell contains a multistage shift register with three All processes inside the cell run in parallel for input data, so that communication between processes is overlapped with computation. Each cell produces its partial result one cycle earlier than the cell to its right . The skew can be accomplished by replacing the register in each cell , which transforms the signal stream, with a multistage shift register. As shown in figure 3, the image data is pumped into the array by the first host, and then send them down to the multifunctional cells. Each multifunctional cell requires 2 input and 2 output channels for communications with the left and right neighbour cells. The data and results are pumped through the array and the accumulated partial products are collected by the final cell , then it sends the final results to the second host. The algorithm is repeated for every input data for each cell concurrently. Only cells on the array boundaries are permitted to communicate with the host and each of the cells communicates with the left and right neighbour cells only. Each cell forms a second level of pipelining, the first level being the global pipelining between array cells, while this additional level of pipelining of computations within a cell can increase the system throughput. In each step, all cells simultaneously perform their I/O and execute their operations.

PVM implementation for the systolic design

The systolic system described previously required further modification in order for the array to be implemented on distributed computing using PVM. The master-slave model in which a separate 'control' program termed the master is responsible for process spawning, initialisation, collection and display of results. The slave programs perform the actual computation (cell) involved; they either are allocated their workloads by the master or perform the allocations themselves. The first and the final slaves are connected to the master (host), the first slave receives all the input data from the master and pump it to the network, the final slave receives the results from the neighbouring slave and sends it down to the master. For simplicity each slave is responsible for a cell of the systolic.

Experimentations

For trying, our new applications, we, first, execute them on a network of stations; but due to poor efficiency of Ethernet communications, the results were not good. The next experimentations have been done on a parallel machine (Cray T3d) which implements PVM upon a tore topology [START_REF] Vincent | Applications sur cray t3d. tome 1 : Généralités mpp et t3d[END_REF][START_REF] Vincent | Applications sur cray t3d. : Tome iv : Documentation shmem et pvm 3.3[END_REF] which proposes an architectue with fast links between processors (with 1 s of latency and a bandwidth of 300 Mbytes/sec). We have also improved the performances of systolic implementation by mixing macro-pipeline methods and multi-pipeline implementation.

Grouping messages : Macro-pipeline

The amount of computing is so small in each cell that the application spends most of the execution time to communicate data and to wait pixels for his neighbours. So we have decided to increase the application granularity by using macro-pipeline methods. Instead of sending a message for each pixel, each cell constructs and sends packets of pixels. In table 1, we can see that we have greatly improve the performances of our systolic Laplacian filter by using macro-pipeline methods. For example, on a picture of 128x128 pixels, we can apply our Laplacian filter in 0.94 second with packet of 10 pixels instead of spending 53 seconds with small packets of 2 pixels.

The figure 4 shows more precisely the execution time of the filter, depending on packet size, with a given image size. We can see that the optimal paket size for applying a Laplacian filter on a 32x32 picture is closed to 22 pixels.

Grouping cells : Multi-pipeline

Our main purpose is to developp a complete image processing library on a distributed architecture. Instead of using the same pipeline to apply all the filters, we have implemented a multi-pipeline application. By doing different low-level image processing in parallel on different pipelines, we compute in parallel different pictures. The results presented in the final paper show that the execution time for computing one picture is very closed of the one for obtaining different pictures with multi-pipeline version. The small differences result of initialisingdata and multi-results collection on master slave. Moreover we also demonstrate that merging various filters on the same pipeline of cells can improve the performances of the filtering.

Conclusion and future works

The PVM implementation of pipelined filters have shown interesting results on Cray T3D. But to improve the performances, we have to developp macro-pipeline methods to reduce the communication time of the application. Moreover, the use of multi-pipeline low-level image processing improves greatly the speedup of the application by outputing different pictures at the same time. We are currently finishing the implementation of the whole image processing library with various kinds of operators (convolution, Gamma filter, Inverse Gradient filter...). To compare our systolic view of the problem filtering with other techniques, we are developping also different approachs of the problem by splitting the picture and sending the parts to the processors (data-parallelism approach). We also study a distributed-shared memory implementation to show that these filtering application can be implemented in an easy way on that kind of system.

Figure 3 :

 3 Figure3: Systolic Array for Laplacian different numbers of shift registers in three different multifunctional cells. All processes inside the cell run in parallel for input data, so that communication between processes is overlapped with computation. Each cell produces its partial result one cycle earlier than the cell to its right . The skew can be accomplished by replacing the register in each cell , which transforms the signal stream, with a multistage shift register. As shown in figure3, the image data is pumped into the array by the first host, and then send them down to the multifunctional cells. Each multifunctional cell requires 2 input and 2 output channels for communications with the left and right neighbour cells. The data and results are pumped through the array and the accumulated partial products are collected by the final cell , then it sends the final results to the second host. The algorithm is repeated for every input data for each cell concurrently. Only cells on the array boundaries are permitted to communicate with the host and each of the cells communicates with the left and right neighbour cells

Figure 4 :

 4 Figure4: Laplacian filter on a small 32x32 picture, depending on packet size (pipeline width)

Table 1 :

 1 Laplacian filter execution time depending on packet size for various picture size (in seconds)

		16x16 32x32 64x64 128x128
	1	0.066	0.62	12.1
	2	0.035	0.19	2.5	52.3
	10	0.016	0.06	0.23	0.94
	100 0.048	0.1	0.25	1.01

This work has been supported by The Human Capital & Mobility EC Project NATHAN