/
blaschkemet.py
122 lines (98 loc) · 4.14 KB
/
blaschkemet.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
import numpy as np
import scipy.interpolate
import scipy.integrate
import disclap
from wanglin import WangLinearization
import functools
# CONVENTIONS
# Infinitesimal invariants of the affine sphere are passed around as "a-vectors",
# which are 5-tuples consisting of
# (u(p), du/dx(p), du/dy(p), Re(C(p)), Im(C(p)))
# The affine frame field is considered as a 3x3 matrix in which the ROWS are
# f, f_x, f_y
def coefmat_dx(a):
u = a[0]
ux = a[1]
uy = a[2]
cr = a[3]
ci = a[4]
emu = np.exp(-u)
return np.matrix([ [0.0, 1.0, 0.0],
[np.exp(u), cr*emu + 0.5*ux, -ci*emu - 0.5*uy],
[0.0, -ci*emu + 0.5*uy, -cr*emu + 0.5*ux] ])
def coefmat_dy(a):
u = a[0]
ux = a[1]
uy = a[2]
cr = a[3]
ci = a[4]
emu = np.exp(-u)
return np.matrix([ [0.0, 0.0, 1.0],
[0.0, -ci*emu + 0.5*uy, -cr*emu + 0.5*ux],
[np.exp(u), -cr*emu - 0.5*ux, ci*emu + 0.5*uy] ])
def coefmat(theta,a):
return np.cos(theta)*coefmat_dx(a) + np.sin(theta)*coefmat_dy(a)
def odeA(theta,pcdata,t,y):
m = coefmat(theta,pcdata(t * np.exp(1j*theta)))
ya = np.matrix(y)
ya.shape = (3,3)
return (m * ya).ravel() # flatten 3x3 matrix to 9-vector
def phase(z):
a = np.angle(z)
if a < 0:
return 2.0*np.pi+a
else:
return a
class BlaschkeMetric(WangLinearization):
DEFAULT_RET_TYPE='affine'
def __init__(self,*args,**kwargs):
WangLinearization.__init__(self,*args,**kwargs)
# Prepare for ODE integration
self.ux = disclap.discdx(self.grid) * self.u
self.uy = disclap.discdy(self.grid) * self.u
self.yinit = np.eye(3).ravel()
self.uinterp = scipy.interpolate.RectBivariateSpline(self.grid.x,self.grid.y,np.transpose(self.u.reshape((self.grid.nx,self.grid.ny))))
self.uxinterp = scipy.interpolate.RectBivariateSpline(self.grid.x,self.grid.y,np.transpose(self.ux.reshape((self.grid.nx,self.grid.ny))))
self.uyinterp = scipy.interpolate.RectBivariateSpline(self.grid.x,self.grid.y,np.transpose(self.uy.reshape((self.grid.nx,self.grid.ny))))
def pcdata(self,z):
cval = self.c(z)
return np.array( [self.uinterp(z.real,z.imag),
self.uxinterp(z.real,z.imag),
self.uyinterp(z.real,z.imag),
cval.real,
cval.imag] )
def prepare_output(self,F,return_type='affine'):
if return_type == 'affine':
return F[0,1:] / F[0,0]
elif return_type == 'homog':
return F[0,:]
elif return_type == 'frame':
return F
def integrate_polar(self,r=None,theta=0.0,step=0.01,tol=0.00001,return_type='affine'):
# If r not given, use circle 90% of radius of grid incircle
if r==None:
r = 0.9 * self.grid.r
odef = functools.partial(odeA,theta,self.pcdata)
solver = scipy.integrate.ode(odef)
solver.set_integrator("vode",method="adams",with_jacobian=False,first_step=(step*r),rtol=tol,nsteps=50000)
solver.set_initial_value(self.yinit,0.0)
solver.integrate(r)
return self.prepare_output(solver.y.reshape((3,3)), return_type)
def integrate_point(self,z,step=0.01,tol=0.00001,return_type=DEFAULT_RET_TYPE):
return self.integrate_polar(r=abs(z),theta=np.angle(z),step=step,tol=tol,return_type=return_type)
def integrate_rays(self,n,r=None,theta0=0.0,step=0.01,tol=0.0001,return_type=DEFAULT_RET_TYPE):
thetalist = np.linspace(theta0, theta0 + 2*np.pi, num=n, endpoint=False)
return [ self.integrate_polar(r=r,theta=theta,step=step,tol=tol,return_type=return_type) for theta in thetalist ]
def integrate_points(self,zlist,step=0.01,tol=0.00001,return_type=DEFAULT_RET_TYPE):
return [ self.integrate_point(z,step=step,tol=tol,return_type=return_type) for z in zlist ]
def _moduletest():
import squaregrid
def c(z):
return z*z
gr = squaregrid.SquareGrid(5.0,100)
B = BlaschkeMetric(c,gr)
vlist = B.integrate_rays(50)
for v in vlist:
print v[0],v[1]
if __name__=='__main__':
_moduletest()