diff --git a/README.md b/README.md
index 8029d46..b5fdb31 100644
--- a/README.md
+++ b/README.md
@@ -26,19 +26,34 @@ The problem of *quadratization* is, given a system of ODEs with polynomial right
system with quadratic right-hand side by introducing as few new variables as possible. We will explain it using toy
example. Consider the system
-
+$$
+\begin{cases}
+x_1' = x_1 x_2 \\
+x_2' = -x_1 x_2^3
+\end{cases}
+$$
-An example of quadratization of this system will be a new variable
+An example of quadratization of this system will be a singular vector of new variables
-
+$$
+[y' = x_2 y - 2y^2]
+$$
leading to the following ODE
-
+$$
+y' = x_2 y - 2y^2
+$$
Thus, we attained the system with quadratic right-hand side
-
+$$
+\begin{cases}
+x_1' = x_1 x_2 \\
+x_2 ' = -x^2 y \\
+y' = x_2 y - 2y^2
+\end{cases}
+$$
We used only one new variable, so we achieved an *optimal* quadratization.
@@ -63,7 +78,14 @@ from qbee import *
For example, we will take the **A1** system from [Swischuk et al.'2020](https://arxiv.org/abs/1908.03620)
-
+$$
+\begin{cases}
+c_1' = -A \exp(-E_a / (R_u T)) c_1 ^{0.2} c_2^{1.3} \\
+c_2 ' = 2c_1' \\
+c_3' = -c_1' \\
+c_4' = -2 c_1'
+\end{cases}
+$$
The parameters in the system are `A, Ea and Ru`, and the others are either state variables or inputs.
So, let's define them with the system in code:
diff --git a/docs/build.py b/docs/build.py
new file mode 100644
index 0000000..407694b
--- /dev/null
+++ b/docs/build.py
@@ -0,0 +1,5 @@
+import subprocess
+
+
+def build_html():
+ subprocess.run("sphinx-build -b html docs/source docs/build")
diff --git a/docs/source/articles_and_resources.rst b/docs/source/articles_and_resources.rst
new file mode 100644
index 0000000..c3aa544
--- /dev/null
+++ b/docs/source/articles_and_resources.rst
@@ -0,0 +1,13 @@
+Articles and Resources
+=======================
+
+* **2023** -- `Exact and optimal quadratization of nonlinear finite-dimensional non-autonomous dynamical systems `_
+* **2021** -- `Optimal monomial quadratization for ODE systems `_ (`arxiv `_)
+ * `Talk at IWOCA 2021 `_, `slides `_
+* **2020** -- `Quadratization of ODEs: Monomial vs. Non-Monomial `_
+ * `Talk at Young Mathematicians Conference 2020 `_
+* **2020** -- `Optimal monomial quadratization for ODE systems: extended abstract `_
+ * `Talk at ISSAC 2020 `_, `slides `_
+
+
+#TODO: Add articles of Lifeware, Boris' group, etc.
\ No newline at end of file
diff --git a/docs/source/background.rst b/docs/source/background.rst
new file mode 100644
index 0000000..357c418
--- /dev/null
+++ b/docs/source/background.rst
@@ -0,0 +1,3 @@
+Background
+===================
+
diff --git a/docs/source/cheatsheet.rst b/docs/source/cheatsheet.rst
new file mode 100644
index 0000000..3f60efd
--- /dev/null
+++ b/docs/source/cheatsheet.rst
@@ -0,0 +1,3 @@
+Cheatsheet
+===================
+
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 486ede3..d480098 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -18,7 +18,7 @@
# -- Project information -----------------------------------------------------
project = 'QBee'
-copyright = '2021, Andrey Bychkov, Gleb Pogudin'
+copyright = '2023, Andrey Bychkov, Gleb Pogudin'
author = 'Andrey Bychkov, Gleb Pogudin'
# The full version, including alpha/beta/rc tags
@@ -34,7 +34,9 @@
"sphinx.ext.autodoc",
"sphinx.ext.napoleon",
"sphinx.ext.autosectionlabel",
- "sphinx.ext.mathjax"
+ "sphinx.ext.mathjax",
+ "sphinx_math_dollar",
+ "sphinx_toolbox.collapse"
]
# Add any paths that contain templates here, relative to this directory.
diff --git a/docs/source/dimension_agnostic.rst b/docs/source/dimension_agnostic.rst
new file mode 100644
index 0000000..59a1f8a
--- /dev/null
+++ b/docs/source/dimension_agnostic.rst
@@ -0,0 +1,4 @@
+Dimension-Agnostic Quadratization
+==================================
+
+Lorem ipsum
\ No newline at end of file
diff --git a/docs/source/images/Diff_operators.png b/docs/source/images/Diff_operators.png
new file mode 100644
index 0000000..124ad55
Binary files /dev/null and b/docs/source/images/Diff_operators.png differ
diff --git a/docs/source/index.rst b/docs/source/index.rst
index ecd4128..f5bedc8 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -31,4 +31,11 @@ Indices and tables
:caption: Contents:
installation
+ cheatsheet
+ quadratization
+ inputs_and_parameters
+ polynomialization
+ dimension_agnostic
+ numerical_integration
+ articles_and_resources
api
\ No newline at end of file
diff --git a/docs/source/inputs_and_parameters.rst b/docs/source/inputs_and_parameters.rst
new file mode 100644
index 0000000..d706e47
--- /dev/null
+++ b/docs/source/inputs_and_parameters.rst
@@ -0,0 +1,252 @@
+Inputs and Parameters
+======================
+
+Differential equations may also contain not only state variables, but also parameters and input functions.
+For instance, consider the system:
+$$
+\begin{array}{l}
+x' = c_1 y \cdot u(t) \\
+y' = c_2 x^3
+\end{array},
+$$
+
+where $c_1, c_2 \in \mathbb{C}$ are parameters and $u(t)$ is an unknown input function.
+
+In QBee we introduce them in the following way:
+
+.. code-block:: python
+
+ >>> from qbee import *
+ >>> c1, c2 = parameters("c1, c2")
+ >>> x, y, u = functions("x, y, u")
+ >>> system = [(x, c1*y*u), (y, c2*x**3)] # u is considered as input since it's not in the left-hand side;
+ >>> quadratize(system).print()
+ Introduced variables:
+ w0 = x**2
+ w1 = x*y
+ w2 = y**2
+
+ x' = c1*u*y
+ y' = c2*w0*x
+ w0' = 2*c1*u*w1
+ w1' = c1*u*w2 + c2*w0**2
+ w2' = 2*c2*w0*w1
+
+
+Inputs' Derivatives
+--------------------------
+
+Note that we can introduce a derivative of an input function during quadratization.
+Let's look at this with an example:
+$$
+\begin{array}{ll}
+x' = y^2 u(t) \\
+y' = x^2
+\end{array}
+$$
+
+.. code-block:: python
+
+ >>> from qbee import *
+ >>> x, y, u = functions("x, y, u")
+ >>> system = [(x, y**2 * u), (y, x**2)]
+ >>> quadratize(system).print()
+ Introduced variables:
+ w0 = u*y
+ w1 = u*x
+
+ x' = w0*y
+ y' = x**2
+ w0' = u'*y + w1*x
+ w1' = u'*x + w0**2
+
+It should be mentioned that sometimes the derivatives of the input functions are unknown or do not even exist!
+To cover such cases, we can try to find an *input-free quadratization*,
+i.e. quadratization contains no input functions and therefore does not introduce their derivatives.
+
+.. code-block:: python
+
+ >>> from qbee import *
+ >>> x, y, u = functions("x, y, u")
+ >>> system = [(x, y**2 * u), (y, x**2)]
+ >>> quadratize(system, input_free=True).print()
+ Introduced variables:
+ w0 = y**2
+ w1 = x**2
+ w2 = x*y**2
+ w3 = x*y
+ w4 = y**3
+ w5 = y**4
+
+ x' = u*w0
+ y' = w1
+ w0' = 2*w1*y
+ w1' = 2*u*w2
+ w2' = u*w5 + 2*w1*w3
+ w3' = u*w4 + w1*x
+ w4' = 3*w3**2
+ w5' = 4*w2*w3
+
+It is clear that the order of quadratization in the input-free case has increased significantly.
+This is not accidental: such quadratizations are indeed larger and do not even always exist!
+
+Existence of input-free quadratizations
+---------------------------------------------------------
+
+The existence from problem of input-free quadratization is thoroughly discussed in Section 3.2 of
+`Exact and optimal quadratization of nonlinear finite-dimensional non-autonomous dynamical systems `_ .
+Here, however, we will only give our results and helpful practices.
+
+Consider original system to be input-affine, that is, of the form
+$$
+\mathbf{\dot{x}} = \mathbf{p_0}(\mathbf{x}) + \sum_{i=1}^r \mathbf{p_i}(\mathbf{x}) u_i. \tag{1}
+$$
+
+**Theorem 1**: if the system (1) is also a polynomial-bilinear, that is, the total degree of $\mathbf{p_i}(\mathbf{x})$
+is at most one for every $1 \le i \le r$, then there is an input-free quadratization.
+
+
+It turns out that the existence of input-free quadratization can be characterized via
+the properties of certain linear differential operators associated with inputs $\mathbf{u}(t)$.
+
+**Definition 1**: We introduce $r$ differential operators for the system (1):
+$$
+D_i := \mathbf{p_i}(\mathbf{x})^T \cdot \frac{\partial}{\partial \mathbf{x}}, \quad 1 \le i \le r,
+$$
+where $\frac{\partial}{\partial \mathbf{x}} = \Big[ \frac{\partial}{\partial x_1},\dots \frac{\partial}{\partial x_N} \Big]^T$
+
+It is easier to imagine with an **example**. Consider the following system:
+$$
+\begin{array}{l}
+\mathbf{\dot{x_1}} = x_1 + x_1 y_1 \\
+\mathbf{\dot{x_2}} = x_1^2 u_1 + x_2 u_2
+\end{array}.
+$$
+
+Then we introduce two differential operators $D_1$ and $D_2$, associated with $u_1$ and $u_2$ accordingly:
+
+$$
+D_1 = [x_1, x_1^2] \cdot [\frac{\partial}{\partial x_1}, \frac{\partial}{\partial x_2}]
+$$
+and
+$$
+D_2 = [0, x_2] \cdot [\frac{\partial}{\partial x_1}, \frac{\partial}{\partial x_2}].
+$$
+
+Below you can see how we build up these operators more visually:
+
+.. image:: images/Diff_operators.png
+ :alt: Differential operators construction
+ :width: 600
+
+**Theorem 2**: Let $\mathcal{A}$ be a subalgebrra generated by $D_1,\dots, D_r$ in the algebra
+$\mathbb{C}[\mathbf{x}, \mathbf{\frac{\partial}{\partial x}}]$ of all polynomial
+differential operators in $\mathbf{x}$.
+Then there is an input-free quadratization of (1) if and only if
+$$
+dim\{ A(x_i)\ | \ A \in \mathcal{A}\} < \infty \quad \textit{for every } 1 \le i \le N.
+$$
+
+Let's use the example above to show how the theorem works.
+
+For $x_1$ we have
+$$
+\begin{gather*}
+D_1(x_1) = x_1, \quad D_2(x_1) = 0 \\ \Downarrow \\
+\{A(x_1)\ | \ A \in \mathcal(A)\} = \text{span}\{x_1\}
+\end{gather*}
+$$
+
+For $x_2$ we have
+$$
+\begin{gather*}
+D_1(x_2) = x_1^2, \quad D_1^2(x_2) = D_1(x_1^2) = 2x_1^2,\quad D_2(x_2) = x_2\\
+\Downarrow \\
+\{A(x_2)\ | \ A \in \mathcal(A)\} = \text{span}\{x_2, x_1^2\}
+\end{gather*}
+$$
+
+Therefore, Theorem 2 implies that the example system has an input-free quadratization.
+And indeed:
+
+.. code-block:: python
+
+ >>> x1, x2, u1, u2 = functions("x1, x2, u1, u2")
+ >>> system = [(x1, x1 + x1 * u1), (x2, x1**2 * u1 + x2 * u2)]
+ >>> quadratize(system, input_free=True).print()
+ Introduced variables:
+ w0 = x1**2
+
+ x1' = u1*x1 + x1
+ x2' = u1*w0 + u2*x2
+ w0' = 2*u1*w0 + 2*w0
+
+
+Now consider an **example** where there is no input-free quadratization.
+$$
+\dot x = x^2 u
+$$
+
+Associated differential operator is $D = x^2 \frac{\partial}{\partial x}$ and the algebra $\mathcal{A}$
+spanned by
+$$
+D(x) = x^2,\quad D^2(x) = D(x^2) = 2x^3,\quad D^3(x) = 6x^4, \dots
+$$
+has infinite dimension. Thus, there is no input-free quadratization exists.
+
+.. code-block:: python
+
+ >>> x, u = functions("x, u")
+ >>> system = [(x, x**2 * u)]
+ >>> upper_bound = partial(pruning_by_vars_number, nvars=100) # Max quadratization order = 100
+ >>> res = quadratize(system, input_free=True, pruning_functions=[upper_bound, *default_pruning_rules])
+ >>> print(res)
+ None
+
+.. collapse:: Why Theorem 2 is not in QBee yet
+
+ The problem is that a complete algorithm for solving this kind of problem
+ has not yet been found in the general case (See **Remark 3.10** of the article above).
+ For the input-free quadratization problem specifically, there is some progress,
+ but we have not yet been able to fully prove the correctness of the proposed algorithm.
+
+|
+
+Setting up possible derivatives of inputs
+----------------------------------------------------
+
+Sometimes there are cases where input-free quadratizations are not required
+and you know exactly how many derivatives your input functions have.
+
+You can use the ``input_ders_order`` parameter to denote
+the maximum allowed order of derivatives for input functions.
+
+.. code-block:: python
+
+ >>> x, y, u, v = functions("x, y, u, v")
+ >>> system = [(x, y**2 * u), (y, v * x**2)]
+ >>> quadratize(system, input_der_orders={v: 2, u: 0}).print() # v' and v'' are allowed
+ Introduced variables:
+ w0 = y**2
+ w1 = v*x
+ w2 = v*y**2
+ w3 = v*x*y
+ w4 = v'*x
+ w5 = v'*y**2
+ w6 = v*x**2
+ w7 = v*y**3
+ w8 = v*x*y**2
+ w9 = v*y**4
+
+ x' = u*w0
+ y' = w6
+ w0' = 2*w6*y
+ w1' = u*w2 + w4
+ w2' = 2*w1*w3 + w5
+ w3' = u*w7 + w1*w6 + w4*y
+ w4' = u*w5 + v''*x
+ w5' = v''*w0 + 2*w3*w4
+ w6' = 2*u*w8 + w4*x
+ w7' = 3*w3**2 + w5*y
+ w8' = u*w9 + w0*w4 + 2*w3*w6
+ w9' = w0*w5 + 4*w6*w7
\ No newline at end of file
diff --git a/docs/source/numerical_integration.rst b/docs/source/numerical_integration.rst
new file mode 100644
index 0000000..82a636f
--- /dev/null
+++ b/docs/source/numerical_integration.rst
@@ -0,0 +1,58 @@
+Numerical Integration
+======================
+
+You can quickly check the quadratization numerically by using a function ``qbee.experimental.to_odeint``
+that turns the result into a numerical function that scipy understands. Let's look at the in the example:
+
+We define and quadratize a scalar system
+$$
+x' = p \sin(x) + u(t)
+$$
+
+.. code-block:: python
+
+ >>> import sympy as sp
+ >>> from qbee import *
+ >>> x, u = functions("x, u")
+ >>> p = parameters("p")
+ >>> system = [(x, p * sp.sin(x) + u)]
+ >>> res = polynomialize_and_quadratize(system, input_free=True)
+ >>> res.print(str_function=sp.latex)
+ Introduced variables:
+ w_{0} = \sin{\left(x \right)}
+ w_{1} = \cos{\left(x \right)}
+
+ x' = p w_{0} + u
+ w_{0}' = p w_{0} w_{1} + u w_{1}
+ w_{1}' = - p w_{0}^{2} - u w_{0}
+
+Here we took advantage of the ability to both get a quadratic system in LaTeX and render it in the documentation.
+
+$$
+w_{0} = \sin{\left(x \right)}, w_{1} = \cos{\left(x \right)}
+$$
+$$
+\begin{array}{l}
+x' = p w_{0} + u \\
+w_{0}' = p w_{0} w_{1} + u w_{1} \\
+w_{1}' = - p w_{0}^{2} - u w_{0}
+\end{array}
+$$
+
+Now we solve the system numerically for $t \in [0; 100]$, $p = 0.1$, $u(t) = \sin(t)$:
+
+.. code-block:: python
+
+ >>> from qbee.experimental import to_odeint
+ >>> ts = np.linspace(0, 100, 10000)
+ >>> t = INDEPENDENT_VARIABLE
+ >>> my_odeint = to_odeint(system, {"x": 0.1}, {"p": 0.1}, {"u": sp.sin(t)})
+ >>> my_odeint_quad = to_odeint(res, {"x": 0.1}, {"p": 0.1}, {"u": sp.sin(t)})
+
+ >>> num_res = my_odeint(ts, rtol=1e-10)
+ >>> num_res_quad = my_odeint_quad(ts, rtol=1e-10)
+ >>> print("L-infinity distance between numerical solutions of the original ODE and the quadratized one: ",
+ >>> np.linalg.norm(num_res[:, 0] - num_res_quad[:, 0], np.inf))
+ L-infinity distance between numerical solutions of the original ODE and the quadratized one: 1.5569859268538266e-07
+
+
diff --git a/docs/source/polynomialization.rst b/docs/source/polynomialization.rst
new file mode 100644
index 0000000..d9af394
--- /dev/null
+++ b/docs/source/polynomialization.rst
@@ -0,0 +1,99 @@
+Polynomialization
+===================
+
+Polynomialization by example
+------------------------------
+
+Polynomialization is extremely close in spirit to :ref:`quadratization `.
+The difference is that this time we are converting systems with the right-hand side as elementary functions
+to systems with polynomial right-hand sides.
+Let's look at an example:
+$$
+x' = x \sin(a x) + u(t)
+$$
+
+.. code-block:: python
+
+ >>> import sympy as sp
+ >>> from qbee import *
+ >>> x, u = functions("x, u")
+ >>> a = parameters("a")
+ >>> system = [(x, x * sp.sin(a*x) + u)] # u is an input variable
+ >>> polynomialize(system).print()
+ Introduced variables:
+ w_0 = sin(a*x)
+ w_1 = cos(a*x)
+
+ x' = u + w_0*x
+ w_0' = a*w_1*(u + w_0*x)
+ w_1' = -a*w_0*(u + w_0*x)
+
+Note that the right-hand side of the resulting system is polynomial.
+We will call the set of introduced variables $\{\cos(a \cdot x), \sin(a \cdot x)\}$ a **polynomialization** of the system.
+
+Result of polynomialization can be natively combined with quadratization:
+
+.. code-block:: python
+
+ >>> poly_res = polynomialize(system)
+ >>> quadratize(poly_res, new_vars_name="p_").print()
+ Introduced variables:
+ p_0 = w_1*x
+ p_1 = w_0*x
+
+ x' = p_0 + u
+ w_0' = -a*p_0*w_1 - a*u*w_1
+ w_1' = a*p_1*w_1 + a*u*w_0
+ p_0' = a*p_0*p_1 + a*p_1*u + p_0*w_1 + u*w_1
+ p_1' = -a*p_0**2 - a*p_0*u + p_1*w_1 + u*w_0
+
+Nevertheless, there is a much simpler way to do it altogether:
+
+.. code-block:: python
+
+ >>> polynomialize_and_quadratize(system).print()
+ Introduced variables:
+ w_0 = sin(a*x)
+ w_1 = cos(a*x)
+ w_2 = w_0*x
+ w_3 = w_1*x
+
+ x' = u + w_2
+ w_0' = a*u*w_1 + a*w_1*w_2
+ w_1' = -a*u*w_0 - a*w_0*w_2
+ w_2' = a*u*w_3 + a*w_2*w_3 + u*w_0 + w_0*w_2
+ w_3' = -a*u*w_2 - a*w_2**2 + u*w_1 + w_1*w_2
+
+
+Polynomialization and Laurent monomials
+---------------------------------------------
+
+Consider the following example:
+
+.. code-block:: python
+
+ >>> import sympy as sp
+ >>> from qbee import *
+ >>> x = functions("x")
+ >>> system = [(x, sp.sin(1 / x))]
+ >>> polynomialize(system).print()
+ Introduced variables:
+ w_0 = sin(1/x)
+ w_1 = cos(1/x)
+
+ x' = w_0
+ w_0' = -w_0*w_1/x**2
+ w_1' = w_0**2/x**2
+
+As you can see, the resulting system is polynomial in :ref:`Laurent sense `, i.e.
+have variables with negative integer powers.
+Although this form comes in handy for more optimal quadratization in the end
+(see :ref:`this section `),
+you may want to get a system with a polynomial right-hand side in the usual sense.
+How to do this is shown in the :ref:`corresponding section `.
+
+Optional: Polynomialization of rational powers
+----------------------------------------------
+
+TODO, use my note
+
diff --git a/docs/source/quadratization.rst b/docs/source/quadratization.rst
new file mode 100644
index 0000000..b959fbf
--- /dev/null
+++ b/docs/source/quadratization.rst
@@ -0,0 +1,161 @@
+Quadratization
+================
+
+Quadratization by example
+-------------------------
+
+To better understand quadratization, let's consider the following example. We have a scalar ordinary differential equation (ODE) given by:
+$$
+x' = x^5
+$$
+Then, introduce a new variable $y = x^4$, leading to an ODE with quadratic right-hand side:
+$$
+x' = xy
+$$
+However, this new equation does not provide any additional information.
+To add up, we can differentiate the equation of y.
+This results in a system of ODEs that is equivalent to the original equation:
+$$
+\begin{array}{ll}
+x' = xy \\
+y' = 4x^3 x' = 4x^8 = 4 y^2
+\end{array}
+$$
+
+Moreover, both equations have a quadratic right-hand side.
+We refer to the set $\{y = x^4\}$ as a **quadratization**.
+A formal definition of quadratization will be provided later.
+
+But now, let's use QBee to replicate the example:
+
+.. code-block:: python
+
+ >>> from qbee import *
+ >>> x = functions("x")
+ >>> quadratize([(x, x**5)]).print() # input follows rule (x, p) => x' = p
+ Introduced variables:
+ w0 = x**4
+
+ x' = w0*x
+ w0' = 4*w0**2
+
+Formal definition
+--------------------
+
+Consider a polynomial system of ODEs
+$$
+\mathbf{\dot x} = \mathbf{p}(\mathbf{x}),
+$$
+where $\mathbf{x} = \mathbf{x}(t) = [x_1(t), \dots, x_N(t)]^T$ and $\mathbf{p}(\mathbf{x}) = [p_1(\mathbf{x}),\dots, p_N(\mathbf{x})]^T$
+with $p_1,\dots, p_N \in \mathbb{C} [\mathbf{x}]$ (are polynomials in $\mathbb{C}$).
+
+Then an $l$-dimensional vector of new variables
+$$
+\mathbf{w} = \mathbf{w}(\mathbf{x}) \quad \in \mathbb{C}[\mathbf{x}]^l
+$$
+is said to be a **quadratization** of the given system if there exist vectors $\mathbf{q}_1(\mathbf{x}, \mathbf{w})$
+and $\mathbf{q}_2(\mathbf{x}, \mathbf{w})$ of dimension $N$ and $l$ respectively, with entries being polynomials
+of total degree at most two (quadratic) such that
+$$
+\dot{\mathbf{x}} = \mathbf{q}_1(\mathbf{x}, \mathbf{w}) \quad and \quad \dot{\mathbf{w}} = \mathbf{q}_2(\mathbf{x}, \mathbf{w}).
+$$
+
+The dimension $l$ of vector $\mathbf{w}$ is called **order of quadratization**.
+A quadratization of the smallest possible order is called an **optimal quadratization**.
+
+In the example above $\mathbf{w} = [x^4],\ l = 1,\ \text{and}\ \mathbf{q}_1 = [w_0 x],\ \mathbf{q}_2 = [4w_0^2]$.
+Since we introduced only one variable, the resulting quadratization is optimal.
+
+
+Monomial quadratization
+------------------------
+
+If all the polynomials $w_1(\mathbf{x}),\dots,w_l(\mathbf{x})$ are
+monomials, the quadratization is called a **monomial quadratization**.
+Accordingly, if a monomial quadratization of a system has the smallest possible order among all
+the monomial quadratizations of the system, it is called an **optimal monomial quadratization**.
+
+The quadratization in the first example is indeed an optimal monomial quadratization.
+
+The importance of monomial quadratizations is that optimal monomial quadratizations
+always exist when there is no :ref:`input functions ` involved
+and $w_1(\mathbf{x}),\dots,w_l(\mathbf{x})$ are non-negative monomials
+(:ref:`Otherwise `, quadratization may have a smaller order
+but we are not aware of algorithm to find an optimal quadratization in this case).
+
+The polynomial quadratization, on the other hand, is generally of a lower order, as seen in the example below.
+However, the general algorithm of finding polynomial quadratizations is unknown to us.
+Nevertheless, shifting variables can help you to produce better order when using QBee.
+
+**Example**: The system
+$$
+x' = (x + 1)^{10}
+$$
+has a polynomial quadratization of order $1$:
+$$
+\mathbf{w} = [(x+1)^9],
+$$
+while its optimal monomial quadratization is
+$$
+\mathbf{v} = [x^2, x^4, x^6, x^8, x^9].
+$$
+Replacing $x$ with $y = x+1$ will have a result similar to the polynomial case, yet reachable by using QBee:
+$$
+y' = y^{10}, \quad \mathbf{v} = [x^9].
+$$
+
+As a further read, we recommend an article by Foyez Allaudin "`Quadratization of ODEs: Monomial vs. Non-Monomial `_".
+
+
+Laurent monomials
+----------------------
+
+Laurent monomials are monomials that may be built with both positive and negative powers, for instance $m(x, y) = x^2 y^{-3}$.
+By default, QBee introduces new variables as Laurent monomials, since it results in more optimal quadratizations than by using ordinary monomials.
+
+**Example**: Let's quadratize the system below by using QBee:
+$$
+x' = \exp(1/x).
+$$
+
+.. code-block:: python
+
+ >>> import sympy as sp
+ >>> from qbee import *
+ >>> x = functions("x")
+ >>> polynomialize_and_quadratize([(x, sp.exp(1/x))]).print()
+ Introduced variables:
+ w_0 = exp(1/x)
+ w_1 = w_0/x**2 <- It's a Laurent monomial;
+ w_2 = w_0/x**3 <- This one as well;
+
+ x' = w_0
+ w_0' = -w_0*w_1
+ w_1' = -2*w_0*w_2 - w_1**2
+ w_2' = -3*w_1**2 - w_1*w_2
+
+**Note**: We first polynomialize the system to bring it to a polynomial form. Read about this technique in :ref:`Polynomialization section `.
+
+However, in some contexts, using Laurent monomials is undesirable or even forbidden.
+For such cases, we can disable them:
+
+.. code-block:: python
+
+ >>> import sympy as sp
+ >>> from qbee import *
+ >>> x = functions("x", laurent=False) # Forbid "x" to be of negative powers in monomials
+ >>> polynomialize_and_quadratize([(x, sp.exp(1/x))]).print()
+ Introduced variables: # Note that the quadratization order has grown;
+ w_0 = exp(1/x)
+ w_1 = 1/x <- We introduce 1/x to tackle with the lack of Laurent monomials;
+ w_2 = w_0*w_1**2
+ w_3 = w_0*w_1**3
+
+ x' = w_0
+ w_0' = -w_0*w_2
+ w_1' = -w_2
+ w_2' = -2*w_0*w_3 - w_2**2
+ w_3' = -3*w_2**2 - w_2*w_3
+
+For more details and context, check Remark 5.1 in
+`Exact and optimal quadratization of nonlinear finite-dimensional non-autonomous dynamical systems `_.
\ No newline at end of file
diff --git a/examples/numerical_integration.py b/examples/numerical_integration.py
index 05de2a5..c6a90d1 100644
--- a/examples/numerical_integration.py
+++ b/examples/numerical_integration.py
@@ -22,3 +22,8 @@
num_res_quad = my_odeint_quad(ts, rtol=1e-10)
print(f"L-infinity distance between numerical solutions of the original ODE and the quadratized one: "
f"{np.linalg.norm(num_res[:, 0] - num_res_quad[:, 0], np.inf)}")
+
+ import matplotlib.pyplot as plt
+ plt.plot(ts, np.abs(num_res - num_res_quad))
+ plt.title("Absolute difference between solutions of original and quadratized systems")
+ plt.savefig("../docs/source/images/num_integration_abs_diff.png")
diff --git a/pyproject.toml b/pyproject.toml
index 407d96c..1d13e1f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
[tool.poetry]
name = "qbee"
-version = "0.8.1"
+version = "0.8.2"
description = "Quadratization of differential equations in python"
authors = ["Andrey Bychkov, Gleb Pogudin"]
license = "MIT"
@@ -34,9 +34,10 @@ pygal = "*"
optional = true
[tool.poetry.group.docs.dependencies]
-sphinx = ">=6.0"
+sphinx = "6.2.1"
sphinx-rtd-theme = "*"
sphinx-math-dollar = "*"
+sphinx_toolbox = "3.4.0"
[build-system]
requires = ["poetry-core"]
@@ -51,4 +52,5 @@ markers = [
[tool.poetry.scripts]
tests = "tests:run_tests"
-benchmarks = "tests:run_benchmarks"
\ No newline at end of file
+benchmarks = "tests:run_benchmarks"
+docs = "docs.build:build_html"
\ No newline at end of file
diff --git a/qbee/polynomialization.py b/qbee/polynomialization.py
index e00ceb5..546c132 100644
--- a/qbee/polynomialization.py
+++ b/qbee/polynomialization.py
@@ -13,7 +13,8 @@
from typing import Iterator, Collection, Iterable
from ordered_set import OrderedSet
from .printer import str_qbee
-from .util import make_derivative_symbol, generate_derivatives, key_from_value, Parameter, progress_bar, monom2PolyElem
+from .util import (make_derivative_symbol, generate_derivatives, key_from_value,
+ Parameter, progress_bar, monom2PolyElem)
class VariablesHolder:
@@ -303,15 +304,12 @@ def __len__(self):
POLY_ALGORITHM_INTERRUPTED = False
-def signal_handler(sig_num, frame):
+def signal_handler_poly(sig_num, frame):
global POLY_ALGORITHM_INTERRUPTED
- print("The algorithm has been interrupted. Returning the current best.")
+ print("The polynomialization algorithm has been interrupted. Returning the current best...")
POLY_ALGORITHM_INTERRUPTED = True
-signal.signal(signal.SIGINT, signal_handler)
-
-
def eq_list_to_eq_system(system: list[(sp.Symbol, sp.Expr)]) -> EquationSystem:
lhs, rhs = zip(*system)
# making constant right-hand sides symbolic expressions
@@ -382,7 +380,9 @@ def polynomialize(system: EquationSystem | list[(sp.Symbol, sp.Expr)], upper_bou
system = eq_list_to_eq_system(system)
system.variables.base_var_name = new_vars_name
system.variables.start_new_vars_with = start_new_vars_with
+ signal.signal(signal.SIGINT, signal_handler_poly) # Add Ctrl+C stopper
nvars, opt_system, traversed = poly_algo_step(system, upper_bound, math.inf)
+ final_iter()
return make_laurent(system, opt_system)
@@ -426,7 +426,8 @@ def poly_algo_step(part_res: EquationSystem, upper_bound, best_nvars) -> (int, E
@progress_bar(is_stop=True)
def final_iter():
- pass
+ global POLY_ALGORITHM_INTERRUPTED
+ POLY_ALGORITHM_INTERRUPTED = False
def next_gen(system: EquationSystem) -> Iterator[EquationSystem]:
diff --git a/qbee/quadratization.py b/qbee/quadratization.py
index f93e893..0fd1277 100644
--- a/qbee/quadratization.py
+++ b/qbee/quadratization.py
@@ -227,6 +227,8 @@ def quadratize_poly(polynomials: list[PolyElement],
w0' = w0**2 + w0*w1
w1' = w0*w1 + 2*w1**2
"""
+ signal.signal(signal.SIGINT, signal_handler_quad)
+
if pruning_functions is None:
pruning_functions = default_pruning_rules
@@ -524,15 +526,12 @@ def _final_iter(self):
QUAD_ALGORITHM_INTERRUPTED = False
-def signal_handler(sig_num, frame):
+def signal_handler_quad(sig_num, frame):
global QUAD_ALGORITHM_INTERRUPTED
- print("The algorithm has been interrupted. Returning the current best.")
+ print("The quadratization algorithm has been interrupted. Returning the current best...")
QUAD_ALGORITHM_INTERRUPTED = True
-signal.signal(signal.SIGINT, signal_handler)
-
-
class BranchAndBound(Algorithm):
def domination_upper_bound(self):
@@ -583,6 +582,9 @@ def next_gen(self, part_res: PolynomialSystem):
def _final_iter(self):
self._nodes_traversed = 0
+ global QUAD_ALGORITHM_INTERRUPTED
+ QUAD_ALGORITHM_INTERRUPTED = False
+
@dump_results
def _save_results(self, opt_system):
return [opt_system, ]
@@ -646,9 +648,11 @@ def pruning_by_vars_number(_: Algorithm, system: PolynomialSystem, *args, nvars:
return True
return False
+
def is_var_bound(f):
return type(f) == partial and f.func == pruning_by_vars_number
+
def pruning_by_best_nvars(a: Algorithm, part_res: PolynomialSystem, *args):
"""Branch-and-Bound default pruning """
best_nvars, *_ = args
diff --git a/qbee/util.py b/qbee/util.py
index b6c94c9..2aa7354 100644
--- a/qbee/util.py
+++ b/qbee/util.py
@@ -93,7 +93,7 @@ def repl(f):
@parametrized
def progress_bar(func, is_stop):
- postfix_str = "Current best order = {}"
+ postfix_str = "Current best order = {}, To return current best press Ctrl+C"
@wraps(func)
def wrapper(*args, **kwargs):