agd.ExportedCode.Notebooks_Algo.SternBrocot

 1# Code automatically exported from notebook SternBrocot.ipynb in directory Notebooks_Algo
 2# Do not modify
 3import numpy as np
 4from matplotlib import pyplot as plt
 5
 6def MakeStencil(refine_pred):
 7    l = [np.array([1,0]),np.array([0,-1]),np.array([-1,0]),np.array([0,1])]
 8    m = [np.array([1,0])]
 9    while len(l)>0:
10        u=m[-1]
11        v=l[-1]
12        if(refine_pred(u,v)):
13            l.append(u+v)
14        else:
15            m.append(v)
16            l.pop()
17    return m
18
19def PlotStencil(stencil):
20    plt.plot(*np.array(stencil).T)
21    plt.scatter(*np.array(stencil).T)
22    plt.scatter(0,0,color='black')
23
24aX0 = np.linspace(-1,1); aX1=aX0
25X = np.array(np.meshgrid(aX0,aX1,indexing='ij'))
26
27def ball_and_stencil(metric,level,name):
28    plt.figure(figsize=[12,4])
29    plt.subplot(1,2,1); plt.title("Unit ball for a norm of "+name+" type"); plt.axis('equal')
30    plt.contourf(*X,metric.norm(X),levels=[0.,level]); plt.scatter(0,0,color='black'); 
31    plt.subplot(1,2,2); plt.title("Stencil for a norm of "+name+" type"); plt.axis('equal')
32    PlotStencil(MakeStencil(lambda u,v: metric.angle(u,v)>np.pi/3))
def MakeStencil(refine_pred):
 7def MakeStencil(refine_pred):
 8    l = [np.array([1,0]),np.array([0,-1]),np.array([-1,0]),np.array([0,1])]
 9    m = [np.array([1,0])]
10    while len(l)>0:
11        u=m[-1]
12        v=l[-1]
13        if(refine_pred(u,v)):
14            l.append(u+v)
15        else:
16            m.append(v)
17            l.pop()
18    return m
def PlotStencil(stencil):
20def PlotStencil(stencil):
21    plt.plot(*np.array(stencil).T)
22    plt.scatter(*np.array(stencil).T)
23    plt.scatter(0,0,color='black')
aX0 = array([-1. , -0.95918367, -0.91836735, -0.87755102, -0.83673469, -0.79591837, -0.75510204, -0.71428571, -0.67346939, -0.63265306, -0.59183673, -0.55102041, -0.51020408, -0.46938776, -0.42857143, -0.3877551 , -0.34693878, -0.30612245, -0.26530612, -0.2244898 , -0.18367347, -0.14285714, -0.10204082, -0.06122449, -0.02040816, 0.02040816, 0.06122449, 0.10204082, 0.14285714, 0.18367347, 0.2244898 , 0.26530612, 0.30612245, 0.34693878, 0.3877551 , 0.42857143, 0.46938776, 0.51020408, 0.55102041, 0.59183673, 0.63265306, 0.67346939, 0.71428571, 0.75510204, 0.79591837, 0.83673469, 0.87755102, 0.91836735, 0.95918367, 1. ])
aX1 = array([-1. , -0.95918367, -0.91836735, -0.87755102, -0.83673469, -0.79591837, -0.75510204, -0.71428571, -0.67346939, -0.63265306, -0.59183673, -0.55102041, -0.51020408, -0.46938776, -0.42857143, -0.3877551 , -0.34693878, -0.30612245, -0.26530612, -0.2244898 , -0.18367347, -0.14285714, -0.10204082, -0.06122449, -0.02040816, 0.02040816, 0.06122449, 0.10204082, 0.14285714, 0.18367347, 0.2244898 , 0.26530612, 0.30612245, 0.34693878, 0.3877551 , 0.42857143, 0.46938776, 0.51020408, 0.55102041, 0.59183673, 0.63265306, 0.67346939, 0.71428571, 0.75510204, 0.79591837, 0.83673469, 0.87755102, 0.91836735, 0.95918367, 1. ])
X = array([[[-1. , -1. , -1. , ..., -1. , -1. , -1. ], [-0.95918367, -0.95918367, -0.95918367, ..., -0.95918367, -0.95918367, -0.95918367], [-0.91836735, -0.91836735, -0.91836735, ..., -0.91836735, -0.91836735, -0.91836735], ..., [ 0.91836735, 0.91836735, 0.91836735, ..., 0.91836735, 0.91836735, 0.91836735], [ 0.95918367, 0.95918367, 0.95918367, ..., 0.95918367, 0.95918367, 0.95918367], [ 1. , 1. , 1. , ..., 1. , 1. , 1. ]], [[-1. , -0.95918367, -0.91836735, ..., 0.91836735, 0.95918367, 1. ], [-1. , -0.95918367, -0.91836735, ..., 0.91836735, 0.95918367, 1. ], [-1. , -0.95918367, -0.91836735, ..., 0.91836735, 0.95918367, 1. ], ..., [-1. , -0.95918367, -0.91836735, ..., 0.91836735, 0.95918367, 1. ], [-1. , -0.95918367, -0.91836735, ..., 0.91836735, 0.95918367, 1. ], [-1. , -0.95918367, -0.91836735, ..., 0.91836735, 0.95918367, 1. ]]])
def ball_and_stencil(metric, level, name):
28def ball_and_stencil(metric,level,name):
29    plt.figure(figsize=[12,4])
30    plt.subplot(1,2,1); plt.title("Unit ball for a norm of "+name+" type"); plt.axis('equal')
31    plt.contourf(*X,metric.norm(X),levels=[0.,level]); plt.scatter(0,0,color='black'); 
32    plt.subplot(1,2,2); plt.title("Stencil for a norm of "+name+" type"); plt.axis('equal')
33    PlotStencil(MakeStencil(lambda u,v: metric.angle(u,v)>np.pi/3))