Integration on manifolds exercise in Sympy!
===========================================
Here I compute the scale factor $J(t)$ for the integration over the faces of a "cube mapped" sphere. That is
$$ F = (u,v,1)^t/\sqrt{u^2 + v^2 + 1}$$

Explanation:
------------
Following http://www.owlnet.rice.edu/~fjones/chap11.pdf

Let $M \in R^n$ be a manifold defined by map $F: A \rightarrow M$, where $A \in R^m, m<n$.
If $g$ is a real-valued function defined on $M$, its integral over $M$ is given by
$$ I = \int_M g(x) dx = \int_A g(F(t)) J(t) dt $$
with
$$ J(t) = \sqrt{det \{(DF)^t (DF)\}}$$
and the Jacobian matrix
$$ DF_{ij} = \partial F_i / \partial t_j$$

For probability densities we have
$$ I = \int_M p_X(x) dx = \int_A p_X(F(t)) J(t) dt = \int_A p_T(t) dt $$
Hence
$$ p_X = \frac{p_T}{J} $$


```python
from sympy import *
import sympy
print(sympy.__version__)
```


```python
u, v = Symbol('u'), Symbol('v')
```


```python
n = sqrt(u*u + v*v + 1.)
```


```python
print(n)
```


```python
x = u/n
y = v/n
z = 1/n
```


```python
import numpy as np
```


```python
DF = np.zeros((3,2), np.object)
for i, a in enumerate([x, y, z]):
    for j, b in enumerate([u, v]):
        DF[i,j] = simplify(diff(a, b))
DF = Matrix(DF)
DF
```


```python
DF.T
```


```python
J = sqrt((DF.T*DF).det())
```


```python
J = simplify(J)
```


```python
J
```


```python
J.subs([(u, 0), (v, 0)]).evalf()
```


```python
q = np.linspace(-1., 1., 21)
print(q)
```


```python
# let's check if the integration works by computing the 
# area of a cube mapped face of the unit sphere.
# Of course the area should come out as 1/6 of the area of the full sphere.
du = (q[1]-q[0])
dv = du
def f(s, t):
    return J.subs([(u, q[s]), (v, q[t])]).evalf()
I = 0.
for s in range(len(q)-1):
    for t in range(len(q)-1):
        # 2D trapezoid quadrature.
        vol = (f(s, t) + f(s+1, t) + f(s,t+1) + f(s+1, t+1))*0.25*du*dv
        I += vol
```


```python
import math
print('Relative to exact area:', I/(4.*math.pi)*6.)
```


```python
from sympy.utilities.codegen import codegen
```


```python
x_, y_, z_ = symbols('x_ y_ z_')
[(c_name, c_code), (h_name, c_header)] = codegen(
    [("cube_map_J", J), ('cube_map_F', [Eq(x_,x), Eq(y_,y), Eq(z_,z)])], "C", "test", header=False, empty=False)
```


```python
print(c_code)
```


```python

```

e9678f990db5e98917c7993be52afa775b019e86
ipynb
Jupyter Notebook
misc/IntegrationOnManifoldsSympy.ipynb
DaWelter/NaiveTrace
a904785a0e13c394b2c221bc918cddb41bc8b175
FSFAP
2018-04-25T08:14:14.000Z
2022-01-29T06:19:16.000Z
misc/IntegrationOnManifoldsSympy.ipynb
DaWelter/NaiveTrace
a904785a0e13c394b2c221bc918cddb41bc8b175
FSFAP
misc/IntegrationOnManifoldsSympy.ipynb
DaWelter/NaiveTrace
a904785a0e13c394b2c221bc918cddb41bc8b175
FSFAP
Qwen/Qwen-72B
1. YES
2. YES