From f7702df6df37eaa82cfc120e98852081206df399 Mon Sep 17 00:00:00 2001
From: Ruizhi Yu <rz.yu@foxmail.com>
Date: Fri, 22 Nov 2024 12:16:55 +0800
Subject: [PATCH] docs: change the advanced usage.

---
 docs/src/advanced.md | 235 +++++++++++++++++++++++++++----------------
 1 file changed, 147 insertions(+), 88 deletions(-)

diff --git a/docs/src/advanced.md b/docs/src/advanced.md
index 6f9746e..ee629bd 100644
--- a/docs/src/advanced.md
+++ b/docs/src/advanced.md
@@ -4,130 +4,189 @@
 
 ## Writing Custom Functions
 Sometimes one may have functions that go beyond the Solverz built-in library. This guide will describe how to create 
-such custom functions in Solverz, so that the functions can be incorporated into numerical simulations. The philosophy 
-of function customization comes from Sympy, it helps to learn the [Sympy basics](https://docs.sympy.org/latest/index.html) 
+such custom functions using Solverz and inform Solverz of their paths, so that the functions can be incorporated into 
+numerical simulations. 
+
+```{note}
+Alternatively, one can directly contribute to the [SolMuseum](https://solmuseum.solverz.org/stable/) 
+library so that 1) others can utilize your models/functions and 2) one avoids configuring the module paths.
+```
+
+```{note}
+The philosophy of function customization comes from Sympy, it helps to learn the [Sympy basics](https://docs.sympy.org/latest/index.html) 
 and read the [Sympy tutorial of custom functions](https://docs.sympy.org/latest/guides/custom-functions.html) for an overview.
+```
+
+```{note}
+In Solverz, the numerical computations are mainly dependent on the prevailing numerical libraries such as numpy and scipy. 
+It is recommended that one first gets familiar with the [numpy](https://numpy.org/doc/stable/user/index.html) and 
+[scipy](https://docs.scipy.org/doc/scipy/tutorial/index.html).
+```
 
-As a motivating example for this document, let's create a custom function class representing the $\min$ function. We 
-want use $\min$ to determine the smaller one of two operands, which can be defined by
+### An Illustrative Example
+
+As a motivating example for this document, let's create a custom function class representing the $\min$ function.
+The $\min$ function is typical in controllers of many industrial applications, which can be defined by
 
 ```{math}
 \min(x,y)=
 \begin{cases}
 x&x\leq y\\
 y& x>y
-\end{cases}
+\end{cases}.
 ```
 
-We also want to extend the function to vector input, that is, capable of finding the element-wise minimum of two vectors.
-To summarize, we shall implement $\min$ that
+To incorporate $\min$ in our simulation modelling, its symbolic and numerical implementations shall be defined. 
+Specifically,
 
-1. evaluates $\min(x,y)$ correctly
-2. can be derived proper derivatives with respect to $x$ and $y$.
+1. a symbolic function `min` can be called to represent the $\min$ function;
+2. the symbolic derivatives of `min` are automatically derived for the Jacobian block parser;
+3. the numerical interfaces are defined so that the `min` function and its derivatives can be correctly evaluated.
 
-However, it is difficult to devise the analytical derivatives of $\min$. We should perform the trick that rewrites 
-$\min(x,y)$ as
+First, we define the numerical interfaces. The derivatives of $\min$ function are
 
 ```{math}
-x*\operatorname{lessthan}(x,y)+y*(1-\operatorname{lessthan}(x,y)).
+\pdv{\min(x,y)}{x}=
+\begin{cases}
+1&x\leq y\\
+0& x>y
+\end{cases},\quad
+\pdv{\min(x,y)}{y}=
+\begin{cases}
+0&x\leq y\\
+1& x>y
+\end{cases}.
 ```
 
-where the $\operatorname{lessthan}(x,y)$ function mathematically denotes the $\leq$ operator and returns 1 if 
-$x\leq y$ else 0. Since $\operatorname{lessthan}(x,y)$ can only be either 1 or 0, the above transformation holds. 
+Let us create a `myfunc.py` file and put the numerical codes there, which will look like
 
-If the derivatives of $\operatorname{lessthan}(x,y)$ with respect to any argument are zero, then we have
-```{math}
-\frac{\partial}{\partial x}\min(x,y)=
-\operatorname{lessthan}(x,y)
+```python
+# myfunc.py
+import numpy as np
+from numba import njit
+
+
+@njit(cache=True)
+def Min(x, y):
+    x = np.asarray(x).reshape((-1,))
+    y = np.asarray(y).reshape((-1,))
+
+    z = np.zeros_like(x)
+
+    for i in range(len(x)):
+        if x[i] <= y[i]:
+            z[i] = x[i]
+        else:
+            z[i] = y[i]
+    return z
+
+
+@njit(cache=True)
+def dMindx(x, y):
+    x = np.asarray(x).reshape((-1,))
+    y = np.asarray(y).reshape((-1,))
+
+    z = np.zeros_like(x)
+
+    for i in range(len(x)):
+        if x[i] <= y[i]:
+            z[i] = 1
+        else:
+            z[i] = 0
+    return z
+
+
+@njit(cache=True)
+def dMindy(x, y):
+    return 1-dMindx(x, y)
 ```
-and
-```{math}
-\frac{\partial}{\partial y}\min(x,y)=
-1-\operatorname{lessthan}(x,y).
+
+where we use `Min` to avoid conflicting with the built-in `min` function, 
+the `@njit(cache)` decorator is used to perform the jit-compilation and hence speed up the numerical codes.
+
+Then let us add the path of `myfunc.py` to Solverz, with, for example, 
+
+```python
+from Solverz import add_my_module
+add_my_module([r'D:\myfunc.py'])
 ```
-Hence, it suffices to have a custom $\operatorname{lessthan}(x,y)$ function that
 
-1. evaluates $\operatorname{lessthan}(x,y)$ correctly
-2. has zero-derivative with respect to $x$ or $y$.
+Then Solverz would load your functions defined in the `myfunc` module. One can also reset the paths by calling 
+`reset_my_module_paths`.
+
+As for the symbolic implementation, let us start by creating a `Min.py` file and subclassing `MulVarFunc` there with
 
-Let us start by subclassing `MulVarFunc`
 ```python
-from Solverz.sym_algebra.functions import MulVarFunc
+from Solverz import MulVarFunc
+
 class Min(MulVarFunc):
     pass
-class LessThan(MulVarFunc):
+
+class dMindx(MulVarFunc):
+    pass
+
+class dMindy(MulVarFunc):
     pass
 ```
-The `MulVarFunc` is the base class of multi-variate functions in Solverz. 
-At this point, `Min` has no behaviors defined on it. To automatically evaluate the `Min` function, we ought to define 
-the **_class method_** `eval()`. `eval()` should take the arguments of the function and return the value 
-$x*\operatorname{lessthan}(x,y)+y*(1-\operatorname{lessthan}(x,y))$:
+
+The `MulVarFunc` is the base class of symbolic multi-variate functions in Solverz. 
+
+At this point, `Min` and its derivatives have no behaviors defined on it. To instruct Solverz in the differentiation 
+rule of `Min` and the numerical implementations, we shall add following lines
 ```python
 class Min(MulVarFunc):
-    @classmethod
-    def eval(cls, x, y):
-        return x * LessThan(x, y) + y * (1 - LessThan(x, y))
-```
-```python
->>> from Solverz import Var
->>> Min(Var('x',0),Var('y',0))
-... x*((x)<=(y)) + y*(1 - ((x)<=(y)))
-```
-To define the differentiation of  `LessThan()`, we have
-```python
-from sympy import Integer
-class LessThan(MulVarFunc):
-    """
-    Represents < operator
-    """
+    arglength = 2
 
-    def _eval_derivative(self, s):
-        return Integer(0)
+    def fdiff(self, argindex=1):
+        if argindex == 1:
+            return dMindx(*self.args)
+        elif argindex == 2:
+            return dMindy(*self.args)
+
+    def _numpycode(self, printer, **kwargs):
+        return (f'myfunc.Min' + r'(' +
+                ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')')
 
-    def _sympystr(self, printer, **kwargs):
-        return '(({op1})<=({op2}))'.format(op1=printer._print(self.args[0]),
-                                           op2=printer._print(self.args[1]))
+
+class dMindx(MulVarFunc):
+    arglength = 2
 
     def _numpycode(self, printer, **kwargs):
-        return r'SolLessThan(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
+        return (f'myfunc.dMindx' + r'(' +
+                ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')')
 
-    def _lambdacode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
 
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
-```
-Here, the `_sympystr()` method defines its string representation:
-```python
->>> LessThan(Var('x'), Var('y'))
-... ((x)<=(y))
-```
-The `_eval_derivative()` method forces the derivatives of `LessThan()` to be zero:
-```python
-from Solverz import iVar
->>> Min(Var('x',0),Var('y',0)).diff(iVar('x'))
-... ((x)<=(y))
-```
-where `iVar` is the internal variable type of Solverz, `diff()` is the method to derive derivatives.
+class dMindy(MulVarFunc):
+    arglength = 2
+
+    def _numpycode(self, printer, **kwargs):
+        return (f'myfunc.dMindy' + r'(' +
+                ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')')
 
-The `_numpycode()` function defines what should `LessThan()` be printed to  in numerical codes. Here, we define the 
-`SolLessThan()` as the numerical implementation of  `LessThan()`. Given array `[0,2,-1]` and `[1,2,3]`:
-```python
->>> import numpy as np
->>> SolLessThan(np.array([0, 2, -1]), np.array([1,2,3]))
-... array([1, 0, 1])
-```
-```{note}
-In Solverz, the numerical computations are mainly dependent on the prevailing numerical libraries such as numpy and scipy. 
-It is recommended that one first gets familiar with the [numpy](https://numpy.org/doc/stable/user/index.html) and 
-[scipy](https://docs.scipy.org/doc/scipy/tutorial/index.html).
 ```
-The implementation of `SolLessThan()` should be put in the `Solverz.num_api.custom_function` module:
+
+where the `fdiff` function should return the derivative of the function, without considering the chain rule, 
+with respect to the `argindex`-th variable; the `_numpycode` functions define the numerical implementations of the 
+functions. Since `myfunc` module has been added to the Solverz path, the numerical implementations can be called by
+`myfunc.Min`.
+
+After finish the above procedures, we can finally use the `Min` function in our simulation modelling, with, for example,
+
 ```python
-@njit(cache=True)
-def SolLessThan(x, y):
-    x = np.asarray(x).reshape((-1,))
-    return (x < y).astype(np.int32)
+from Solverz import Model, Var, Eqn, made_numerical, MulVarFunc
+from Min import Min
+
+m = Model()
+m.x = Var('x', [1, 2])
+m.y = Var('y', [3, 4])
+m.f = Eqn('f', Min(m.x, m.y))
+sae, y0 = m.create_instance()
+ae = made_numerical(sae, y0, sparse=True)
 ```
-The `implements_nfunc()` cannot be omitted and the `njit()` decorator enables the numba-based dynamic compilation for efficiency.
 
+We will have the output
+
+```shell
+>>> ae.F(y0)
+array([1.0, 2.0])
+```