Skip to content

Commit

Permalink
docs: change the advanced usage.
Browse files Browse the repository at this point in the history
  • Loading branch information
@rzyu45
rzyu45 committed Nov 22, 2024
1 parent 50f28aa commit f7702df
Showing 1 changed file with 147 additions and 88 deletions.
235 changes: 147 additions & 88 deletions docs/src/advanced.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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])
```

0 comments on commit f7702df

Please sign in to comment.