-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathDemo_ApproxLF.py
134 lines (104 loc) · 3.63 KB
/
Demo_ApproxLF.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
123
124
125
126
127
128
129
130
131
132
133
134
import time
import numpy as np
import matplotlib.pyplot as plt
from VISolver.Domains.ApproxLF import ApproxLF
from VISolver.Solvers.Euler import Euler
from VISolver.Solvers.HeunEuler import HeunEuler
from VISolver.Solver import Solve
from VISolver.Options import (
DescentOptions, Miscellaneous, Reporting, Termination, Initialization)
from VISolver.Log import PrintSimResults, PrintSimStats
from IPython import embed
def Demo():
# __APPROXIMATE_LINEAR_FIELD__##############################################
# Load Dummy Data
X = np.random.rand(1000,2)*2-1
scale = np.array([0.8,1.2])
y = 0.5*np.sum(scale*X**2,axis=1)
# Construct Field
ALF = ApproxLF(X=X,dy=y,eps=1e-8)
# Set Method
# Method = Euler(Domain=ALF,FixStep=True)
Method = HeunEuler(Domain=ALF,Delta0=1e-2)
# Initialize Starting Field
A = np.array([[0,1],[-1,0]])
# A = np.eye(LF.XDim)
# A = np.random.rand(LF.XDim,LF.XDim)
# A = np.array([[5,0.],[0.,5]])
b = np.zeros(ALF.XDim)
Start = np.hstack([A.flatten(),b])
print(Start)
# Set Options
Init = Initialization(Step=-1.0)
Term = Termination(MaxIter=1000) #,Tols=[(LF.error,1e-10)])
Repo = Reporting(Requests=[ALF.error, 'Step', 'F Evaluations',
'Data'])
Misc = Miscellaneous()
Options = DescentOptions(Init,Term,Repo,Misc)
# Print Stats
PrintSimStats(ALF,Method,Options)
# Start Solver
tic = time.time()
Results = Solve(Start,Method,ALF,Options)
toc = time.time() - tic
# Print Results
PrintSimResults(Options,Results,Method,toc)
error = np.asarray(Results.PermStorage[ALF.error])
params = Results.PermStorage['Data'][-1]
A,b = ALF.UnpackFieldParams(params)
print(A)
print(b)
print(error[-1])
# # Initialize Field
# # A = 10*np.random.rand(LF.XDim,LF.XDim)
# # A = (A+A.T)/2
# A = np.array([[0,1],[-1,0]])
# # A = np.eye(LF.XDim)
# # b = 10*np.random.rand(LF.XDim)
# b = np.zeros(LF.XDim)
# # b = np.ones(2)
# # Compute Path Integral
# t = 1
# t0 = 0
# tf = 1
# x0 = np.zeros(LF.XDim)
# # xf = np.ones(LF.XDim)
# xf = np.zeros(LF.XDim)
# y0 = 0
# y1 = LF.predict([A,b],t0,x0,y0,tf,xf,t=t/2)
# print('y(xt(t/2))=',y1)
# y1 = LF.predict([A,b],t0,x0,y0,tf,xf,t=t)
# print('y(xf)=',y1)
# xt = LF.x(np.linspace(0.0001,1,50),t0,tf,x0,xf,[A,b])
# plt.plot(x0[0],x0[1],'*')
# plt.plot(xt[:,0],xt[:,1],'o')
# plt.plot(xf[0],xf[1],marker=(5,1,0),markersize=20)
# plt.xlim([-1,2])
# plt.ylim([-1,2])
# L = LF.Lagrangian(t,t0,x0,y0,tf,xf,[A,b])
# print('Lagrangian=',L)
# EL = LF.EulerLagrange(np.linspace(0.0001,1,10),t0,x0,y0,tf,xf,[A,b])
# print('EL==0?: ',np.allclose(EL,0))
# S = LF.Action(tf,t0,x0,y0,tf,xf,[A,b])
# print('Action=',S)
# fdG = LF.findiff([A,b],t0,x0,y0,tf,xf,t=t)
# print(fdG)
# print('\n')
# G = LF.gradient([A,b],t0,x0,y0,tf,xf,t=t)
# print(G)
# X, Y = np.meshgrid(np.arange(-1, 2, .1), np.arange(-1, 2, .1))
# # points = np.asarray(list(zip(X.ravel(),Y.ravel())))
# points = np.vstack((X.ravel(),Y.ravel())).T
# vectors = LF.Field([A,b],points)
# U = vectors[:,0].reshape(X.shape)
# V = vectors[:,1].reshape(Y.shape)
# Q = plt.quiver(X[::3, ::3], Y[::3, ::3], U[::3, ::3], V[::3, ::3],
# pivot='mid', color='r', units='inches')
# qk = plt.quiverkey(Q, 0.5, 0.03, 1, r'$1 \frac{m}{s}$',
# fontproperties={'weight': 'bold'})
# plt.plot(X[::3, ::3], Y[::3, ::3], 'k.')
# plt.show()
# # print(xt)
# embed()
if __name__ == '__main__':
Demo()